text:
import inspect
import os

import IPython
import numpy as np

def get_modified_times(dirs, f):
    """Return the times at which the `file` in the directories `dirs` was modified.

    Args:
        dirs (str): Path
        file (str): File name

    Returns:
        list: List of modification times
    """
    mtimes = []
    for d in dirs:
        path = os.path.join(d, f)
        if os.path.exists(path):
            mtimes.append(os.path.getmtime(path))
        else:
            mtimes.append(0)
    return mtimes

def get_parent(directory, parent_folder_name):
    """Get the absolute path of some parenting folder of a directory.

    Args:
        directory (str): The path to search.
        parent_folder_name (str): Name of the parenting folder.

    Returns:
        str, None: Directory of parenting folder. None of not found on path.
    """
    if parent_folder_name not in directory.split(os.sep):
        return None
    parent_directory = ['/']
    for d in directory.split(os.sep):
        if d != parent_folder_name:
            parent_directory.append(d)
        if d == parent_folder_name:
            parent_directory.append(d)
            break
    return os.path.join(*parent_directory)

def longest_common_suffix(list_of_paths):
    reversed_paths = [os.path.join(*s.split(os.sep)[::-1]) for s in list_of_paths]
    reversed_lcs = os.path.commonprefix(reversed_paths)
    lcs = os.path.join(*reversed_lcs.split(os.sep)[::-1])
    if lcs:
        return os.path.join(os.sep, lcs)
    else:
        return None

max_stars_repo_name: JakobHavtorn/es-rl

max_issues_repo_name: JakobHavtorn/es-rl

max_forks_repo_name: JakobHavtorn/es-rl


text:
import pytest

import numpy as np
np.random.seed(1)

from boaf.algorithms.mixture import MixtureModel, GMM
from boaf.base_distributions.multivariate import NIW

@pytest.mark.usefixtures("clustering_data")
def test_mixture(clustering_data):

    nclusters, (X, inds) = clustering_data

    # Options as nested dict allows future JSON options
    opts = {
        'nclusters': nclusters,
        'niters': 200,
        'prior':{
            'alpha': np.ones((nclusters))/nclusters,
            'cluster':{
                'nu':2,
                'kappa':1,
                'mu':np.zeros((1,2)),
                'sigma':np.eye(2)
            }
        }
    }
    model = MixtureModel(opts, base_distribution=NIW)
    model.learn(X)
    inds_pred = model.predict(X)

    assert(np.all(inds_pred < nclusters))
    assert(np.all(inds.shape == inds_pred.shape))

    opts['init'] = 'kpp'
    model = MixtureModel(opts, base_distribution=NIW)
    model.learn(X)
    inds_pred = model.predict(X)

    assert(np.all(inds_pred < nclusters))
    assert(np.all(inds.shape == inds_pred.shape))

    opts['init'] = 'rand'
    model = MixtureModel(opts, base_distribution=NIW)
    model.learn(X)
    inds_pred = model.predict(X)

    assert(np.all(inds_pred < nclusters))
    assert(np.all(inds.shape == inds_pred.shape))

@pytest.mark.usefixtures("clustering_data")
def test_GMM(clustering_data):

    nclusters, (X, inds) = clustering_data

    # Options as nested dict allows future JSON options
    opts = {
        'nclusters': nclusters,
        'niters': 200,
        'prior':{
            'alpha': np.ones((nclusters))/nclusters,
            'cluster':{
                'nu':2,
                'kappa':1,
                'mu':np.zeros((1,2)),
                'sigma':np.eye(2)
            }
        }
    }
    model = GMM(opts)
    model.learn(X)
    inds_pred = model.predict(X)
    print(inds_pred)

    assert(np.all(inds_pred < nclusters))
    assert(np.all(inds.shape == inds_pred.shape))

ext: py

max_stars_repo_name: TimothyRogers/BOAF

max_issues_repo_name: TimothyRogers/BOAF

max_forks_repo_name: TimothyRogers/BOAF

reason: import numpy


text:
import os
import datetime as datetime
from netCDF4 import Dataset
import numpy as np
import glob
import string
from matplotlib.pyplot import cm
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
from pylab import *

__author__   = 'Trond Kristiansen'
__email__    = 'me (at) trondkristiansen.com'
__created__  = datetime.datetime(2019, 3, 6)
__modified__ = datetime.datetime(2019, 3, 6)
__version__  = "0.1"
__status__   = "Development"

def getLocations(mygridfile):
	maxPositions=4

	xpos=np.zeros(shape=(maxPositions,1)); ypos=np.zeros(shape=(maxPositions,1)); xynames=np.zeros(shape=(maxPositions,1))
	lonpos=np.zeros(shape=(maxPositions,1)); latpos=np.zeros(shape=(maxPositions,1))
	xpos=[90,130,150,185]
	ypos=[30,25,90,35]
	xynames=["Soerlige Nordsjoen","Tyskebukten","Nordlige Nordsjoen","Skagerrak"]

	cdf = Dataset(mygridfile)
	lats=cdf.variables["lat_rho"][:]
	lons=cdf.variables["lon_rho"][:]
	lonpos = [lons[ypos[i],xpos[i]] for i in xrange(len(xpos[:]))]
	latpos = [lats[ypos[i],xpos[i]] for i in xrange(len(xpos[:]))]

	print "\nSystem setup with station locations:"
	for name,x,y,lo,la in zip(xynames,xpos,ypos,lonpos,latpos):
		print "%s => longitude=%3.2f  latitude=%3.2f (x=%s,y=%s)"%(name,lo,la,x,y)

	return xpos,ypos,xynames,lonpos,latpos

def getIndicesOfLocationsInGLORYS2V3Grid(lonpos,latpos,xynames,myglorys2v3gridfile):

	print "Creating indices for GLORYS2V3 grid (func:getIndicesOfLocationsInGLORYS2V3Grid)"

	cdf = Dataset(myglorys2v3gridfile)
	latsglorys=cdf.variables["nav_lat"][:]
	lonsglorys=cdf.variables["nav_lon"][:]
	numberOfPoints=1
	glorysIndicies=[]
	glorys=np.zeros(shape=(np.shape(lonpos)[0],2))

	for st_lon, st_lat in zip(lonpos,latpos):
		indices, dis = getStationIndices(st_lon, st_lat,lonsglorys,latsglorys,numberOfPoints)
		glorysIndicies.append(indices)

	for i in xrange(len(glorysIndicies)):
		glorys[i,0]=int(glorysIndicies[i][0])
		glorys[i,1]=int(glorysIndicies[i][1])

	return glorysIndicies

def getStationIndices(st_lon,st_lat,longitudes,latitudes,numberOfPoints):

    distance = np.zeros((longitudes.shape),dtype=np.float64)
    listd=[]
    for eta in range(len(latitudes[:,0])):
        for xi in range(len(latitudes[0,:])):
            distance[eta,xi] = np.sqrt( (latitudes[eta,xi]-st_lat)**2.0 + (longitudes[eta, xi] - st_lon)**2.0 )
            listd.append(distance[eta,xi])

    listsIndexes=[]
    listd.sort()
    for i in range(numberOfPoints):
        value=listd[0]
        itemindex=np.where(distance==value)
        listsIndexes.append(itemindex)

        listd.pop(0)

    print ''
    print '=====getStationIndices======'
    print 'Looking for longitude [%3.3f] and latitude [%3.3f]'%(st_lon,st_lat)
    print 'Result ===>'
    for i in range(numberOfPoints):
        print 'Found index pair in gridfile',listsIndexes[i]
        print 'Index corresponds to longitude [%3.3f] and latitude [%3.3f]'%(longitudes[listsIndexes[i][0],listsIndexes[i][1]],latitudes[listsIndexes[i][0],listsIndexes[i][1]])
    print '======================'
    print ''

    dis=[]

    for i in range(numberOfPoints):
        dis.append(np.sqrt( (latitudes[listsIndexes[i][0],listsIndexes[i][1]]-st_lat)**2.0 + (longitudes[listsIndexes[i][0],listsIndexes[i][1]] - st_lon)**2.0 ))

    return listsIndexes[0], dis

def getValuesFromGLORYS(glorysIndicies,myglorys2v3pathT,myglorys2v3pathS,xynames,startYear,debug):

	print "Extracting values from GLORYS2V3 (func:getValuesFromGLORYS)"

	mypaths=[myglorys2v3pathT,myglorys2v3pathS]
	variables=["votemper","vosaline"]
	first=True
	counter=0; varCounter=0

	for variable,datapath in zip(variables,mypaths):

		argument="%s*.nc"%(datapath)
		allFiles = glob.glob(argument)
		allFiles.sort()
		allFilesFiltered=[]

		# Remove files older than startYear (we only want 2009-2012)
		for afile in allFiles:
			(head, filename) = os.path.split(afile)
			l=string.split(filename,'_')
			date = str(l[2])
			if (startYear <= int(date[0:4])):
				allFilesFiltered.append(afile)
		allFiles = []
		allFiles = allFilesFiltered
		allFiles.sort()

		if first:
			allValues=np.zeros(shape=(len(glorysIndicies),len(allFiles),len(variables)))
			allDates=[]
			print "Sorting %s files found in GLORYS2V3 datadirectory"%(len(allFiles))
			first=False

		refdate=datetime.datetime(1948,1,1)

		for afile in allFiles:
			(head, filename) = os.path.split(afile)
			l=string.split(filename,'_')
			date = str(l[2])
			mydate=datetime.datetime(int(date[0:4]),int(date[4:6]),int(date[6:8]))
			if (varCounter==0):
				allDates.append(mydate)

			cdf = Dataset(afile)

			for station in xrange(len(glorysIndicies)):
				data = cdf.variables[variable][0,0,int(glorysIndicies[station][0]),int(glorysIndicies[station][1])]
				if debug:
					print "%s,%i : %s  - data for station (%s,%s) %s = > %s"%(counter,varCounter,mydate,int(glorysIndicies[station][0]),int(glorysIndicies[station][1]),variable,data)
				allValues[station,counter,varCounter]=data
			cdf.close()
			counter+=1
		counter=0
		varCounter+=1
	return allValues, allDates

def getValuesFromNS8KM(xpos,ypos,xynames,datapath,mypattern,debug):

	print "Extracting values from NS8KM (func:getValuesFromNS8KM)"

	first=True
	counter=0
	variables=["temp","salt"]
	argument="%s%s"%(datapath,mypattern)
	allFiles = glob.glob(argument)
	allFiles.sort()

	print "Sorting %s files found in NS8KM datadirectory"%(len(allFiles))

	refdate=datetime.datetime(1948,1,1)

	for afile in allFiles:

		cdf = Dataset(afile)
		times = dataSST = cdf.variables["ocean_time"][:]
		if first:
			allValues=np.zeros(shape=(len(xpos),len(allFiles)*len(times),len(variables)))
			allDates=[]
			first = False

		for t,time in enumerate(times):

			mydate=(refdate + datetime.timedelta(seconds=int(time)))
			allDates.append(mydate)

			for station in xrange(len(xpos)):
				dataSST = cdf.variables["temp"][t,34,int(ypos[station]),int(xpos[station])]
				dataSSS = cdf.variables["salt"][t,34,int(ypos[station]),int(xpos[station])]
				if debug:
					print "%s : temp  - data for station (%s,%s) %s = > %s"%(counter,mydate,int(xpos[station]),int(ypos[station]),dataSST)
					print "%s : salt  - data for station (%s,%s) %s = > %s"%(counter,mydate,int(xpos[station]),int(ypos[station]),dataSSS)
				allValues[station,counter,0]=dataSST
				allValues[station,counter,1]=dataSSS
			counter+=1
		cdf.close()

	return allValues, allDates

def setupSubPlot(subplotIndex):
	from matplotlib.dates import  DateFormatter
	ax=plt.subplot(2, 1, subplotIndex)
	ax.tick_params(axis='both', which='major', labelsize=6)
	ax.tick_params(axis='both', which='minor', labelsize=6)

  	plt.xticks(rotation=45)
  	ax.xaxis.set_major_formatter( DateFormatter('%Y-%m-%d') )
  	return ax

def createTimeseriesPlot(stations,dataNS8KM,datesNS8KM,dataGLORYS2V3,datesGLORYS2V3,startDate,endDate,plotfileName):

	print "Creating timeseries plot (func:createTimeseriesPlot)"
	dates1 = matplotlib.dates.date2num(datesNS8KM)
  	dates2 = matplotlib.dates.date2num(datesGLORYS2V3)

	ax = setupSubPlot(subplotIndex=1)
	fmtNS=["r-o","g-o","b-o","m-o"]
	fmtGL=["r-.","g-.","b-.","m-."]
	plots=[]
  	for i,station in enumerate(stations):
  		# SST
		myplt,=plt.plot_date(dates1, dataNS8KM[i,:,0],fmt=fmtNS[i],markersize=3,xdate=True,ydate=False)
		plots.append(myplt)

	for i,station in enumerate(stations):
  		# SST
  		plt.plot_date(dates2, dataGLORYS2V3[i,:,0],fmt=fmtGL[i],linewidth=2, xdate=True,ydate=False)

  	legend(plots,stations,loc=3,prop={'size':6})
	ax.set_xlim(matplotlib.dates.date2num(startDate),matplotlib.dates.date2num(endDate))

	plt.ylabel('SST')

	ax2 = setupSubPlot(subplotIndex=2)
  	for i,station in enumerate(stations):
  		# SSS
		plt.plot_date(dates1, dataNS8KM[i,:,1],fmt=fmtNS[i],markersize=3,xdate=True,ydate=False)

	for i,station in enumerate(stations):
  		# SSS
		plt.plot_date(dates2, dataGLORYS2V3[i,:,1],fmt=fmtGL[i],linewidth=2, xdate=True,ydate=False)

	plt.ylabel('SSS')
	startLim=matplotlib.dates.date2num(startDate)
	endLim=matplotlib.dates.date2num(endDate)
	print startLim, endLim
	ax2.set_xlim(startLim,endLim)

	if os.path.exists(plotfileName): os.remove(plotfileName)
	plt.savefig(plotfileName,dpi=300)
	print 'Saved figure file %s\n'%(plotfileName)

def main():

	print "Starting program (func:main)"
	debug = False

	gridfile="/Users/trondkr/Dropbox/NIVA/A20/Grid/A20niva_grd_v1.nc"
	xpos,ypos,xynames,lonpos,latpos = getLocations(gridfile)
	glorysIndicies = getIndicesOfLocationsInGLORYS2V3Grid(lonpos,latpos,xynames,myglorys2v3gridfile)
	allGLORYS2V3Values, allGLORYS2V3Dates = getValuesFromGLORYS(glorysIndicies,myglorys2v3pathT,myglorys2v3pathS,xynames,mystartYear,debug)
 	allNS8KMValues, allNS8KMDates = getValuesFromNS8KM(xpos,ypos,xynames,myremotedir,mypattern,debug)

 	# Create plot for date range of GLORYS2V3 defined dates
 	createTimeseriesPlot(xynames,allNS8KMValues,allNS8KMDates,allGLORYS2V3Values,allGLORYS2V3Dates,allGLORYS2V3Dates[0],allGLORYS2V3Dates[-1],"timeseries_NS8KM_vs_GLORYS2V3.png")

 	# Create zoomed plot for date range of NS8KM defined dates
 	createTimeseriesPlot(xynames,allNS8KMValues,allNS8KMDates,allGLORYS2V3Values,allGLORYS2V3Dates,allNS8KMDates[0],allNS8KMDates[-1],"timeseries_NS8KM_vs_GLORYS2V3_zoomed.png")

if __name__ == "__main__":
	main()


text:
from dipsim import util
from dipsim import multiframe
import numpy as np
import matplotlib.pyplot as plt
import os; import time; start = time.time(); print('Running...')

# Main input parameters
row_labels = ['Geometry', r'$\sigma_{\Omega}$', r'Relative Frequency']
nas = [0.7, 0.9, 1.1, 1.3]
col_labels = ['NA='+str(na) for na in nas]

n_pts = 500
n_cols = len(col_labels)
n_rows = 3
inch_fig = 5
dpi = 400

# Generate axes
size = (inch_fig*n_cols, inch_fig*n_rows)
fig, axs = plt.subplots(n_rows, n_cols, figsize=size)
plt.subplots_adjust(wspace=0.2, hspace=0)
if len(col_labels) == 1:
    axs = np.expand_dims(axs, 1)
caxs = util.generate_caxs(axs)

# Compute and plot on axes
for i, na in enumerate(nas):
    exp = multiframe.MultiFrameMicroscope(ill_thetas=[0], det_thetas=[0],
                                          ill_nas=[na], det_nas=[na],
                                          ill_types=['wide'], det_types=['lens'],
                                          colors=['(1,0,0)'], n_frames=4,
                                          n_pts=n_pts, max_photons=1000, n_samp=1.33)

    print('Computing microscope: ' + str(col_labels[i]))
    exp.calc_estimation_stats()

    scene_string = exp.scene_string()
    util.draw_scene(scene_string, my_ax=axs[0,i], dpi=dpi)

    util.plot_sphere(directions=exp.directions, data=exp.sa_uncert,
                     color_norm='log', linthresh=1e-4,
                     color_min=8e-4, color_max=1e0,
                     my_ax=axs[1,i], my_cax=caxs[1,i])

    util.plot_histogram(exp.sa_uncert, ax=axs[2,i])

    caxs[0,i].axis('off')
    caxs[2,i].axis('off')

# Label axes and save
util.label_rows_and_cols(axs, row_labels, col_labels, row_pos=(-0.2, 0.5))
print('Saving final figure.')
fig.savefig('single-crlb.png', dpi=dpi)

print('Total time: '+str(np.round(time.time() - start, 2)))
os.system('say "done"')

max_stars_repo_name: talonchandler/dipsim

max_issues_repo_name: talonchandler/dipsim

max_forks_repo_name: talonchandler/dipsim


text:
import numpy as np

def handle_img_padding(img: np.ndarray, img_zero_crossing: np.ndarray) -> np.ndarray:
    """Runs the Handle Img Padding algorithm

    :param img: The input image. Must be a gray scale image
    :type img: ndarray
    :param img_zero_crossing: The image with zero crossing
    :type img_zero_crossing: ndarray

    :return: Image with padding applied
    :rtype: ndarray
    """

    M1, N1 = img.shape[:2]
    M2, N2 = img_zero_crossing.shape[:2]
    padding_x = int(np.abs(M2 - M1)/2)
    padding_y = int(np.abs(N2 - N1)/2)
    res = img_zero_crossing[padding_x:M1+padding_x, padding_y: N1+padding_y]
    return res


text:
from keras.layers import Dropout
import numpy as np
from keras.layers import TimeDistributed
from keras import backend as K
from keras.layers import Conv2D, MaxPooling2D, Flatten
from keras.layers import Input, LSTM, Dense
from keras.models import Model
from keras.preprocessing.image import ImageDataGenerator
import configs_and_settings
from sklearn.preprocessing import OneHotEncoder
from load_image_data import ImageDataSource
batch_size = 1
# img_width = 4101
img_height = 4101
# img_height = 247
img_width = 247
channels = 1
num_classes = 9

# if K.image_data_format() == 'channels_first':
	# cnn_model_input_tensor_shape = (batch_size, n_channels, image_height, image_width)
	# cnn_model_input_tensor_shape = (frames_per_seq, 3, img_width, img_height)
	# cnn_model_input_tensor_shape = (None, 3, img_width, img_height)
	# cnn_model_input_tensor_shape = (3, img_width, img_height)
	# cnn_model_input_tensor_shape = (frame_channels, img_width, img_height, None)
	# cnn_model_input_tensor_shape = (frame_channels, img_width, img_height, 2)
# 	cnn_model_input_tensor_shape = (frame_channels, img_width, img_height)
# else:
	# (batch_size, image_height, image_width, n_channels)
	# cnn_model_input_tensor_shape = (frames_per_seq, img_width, img_height, 3)
	# cnn_model_input_tensor_shape = (None, img_width, img_height, 3)
	# cnn_model_input_tensor_shape = (img_width, img_height, 3)
	# cnn_model_input_tensor_shape = (None, img_width, img_height, frame_channels)
	# cnn_model_input_tensor_shape = (2, img_width, img_height, frame_channels)
	# cnn_model_input_tensor_shape = (img_width, img_height, frame_channels)

batches_num = 90    # number of batches
frames_num = 164    # number of sequential samples
height = 247
width = 25
channels = 1
input_tensor_shape = (frames_num, height, width, channels)
# cnn_model_input_tensor_shape = (frame_channels, img_width, img_height)
# cnn_model_input_tensor_shape = (img_width, img_height, frame_channels)
# cnn_model_input_tensor_shape = (number_of_frames, img_width, img_height, frame_channels)
# cnn_model_input_tensor_shape = (img_width, img_height, frame_channels)
input_tensor = Input(shape=input_tensor_shape)

conv1 = TimeDistributed(Conv2D(16, (2, 2), activation="relu"))(input_tensor)
# mp1 = TimeDistributed(MaxPooling2D(pool_size=(2, 2)))(conv1)
mp1 = TimeDistributed(MaxPooling2D(pool_size=(2, 2), padding='same'))(conv1)
do1 = TimeDistributed(Dropout(0.5))(mp1)

conv2 = TimeDistributed(Conv2D(32, (2, 2), activation="relu"))(do1)
# mp2 = TimeDistributed(MaxPooling2D(pool_size=(2, 2)))(conv2)
mp2 = TimeDistributed(MaxPooling2D(pool_size=(2, 2), padding='same'))(conv2)
do2 = TimeDistributed(Dropout(0.5))(mp2)

conv3 = TimeDistributed(Conv2D(64, (2, 2), activation="relu"))(do2)
# mp3 = TimeDistributed(MaxPooling2D(pool_size=(2, 2)))(conv3)
mp3 = TimeDistributed(MaxPooling2D(pool_size=(2, 2), padding='same'))(conv3)
# do3 = TimeDistributed(Dropout(0.5))(mp3)

flat = TimeDistributed(Flatten())(mp3)
# flat = TimeDistributed(Flatten())(do3)

lstm = LSTM(256, return_sequences=True, activation='tanh')(flat)
d1 = Dense(128, activation="relu")(lstm)
do4 = Dropout(0.5)(d1)
d2 = Dense(16, activation="relu")(do4)
do5 = Dropout(0.5)(d2)
d3 = Dense(num_classes, activation="softmax")(do5)

model = Model(inputs=input_tensor, outputs=d3)
model.compile(loss='categorical_crossentropy', metrics=['accuracy'], optimizer="adam")

img_data_src = ImageDataSource()
img_dataset = img_data_src.get_right_end_trimmed_imgs_dataset()
X = img_dataset[0]
# X_train = X_train.reshape((X_train.shape[0], 1, X_train.shape[1], X_train.shape[2], X_train.shape[-1]))
y = img_dataset[1]
y = y.reshape((y.shape[0], 1))
# label_encoder = LabelEncoder()
onehot_encoder = OneHotEncoder(sparse=False)
y = onehot_encoder.fit_transform(y)

num_tot_samples = X.shape[0]
train_split = 0.8
# val_split = 0.2
epochs = 10
batch_size = 1

from sklearn.model_selection import train_test_split
val_split = 0.2
X_train, X_val, y_train, y_val = train_test_split(X, y, test_size=val_split, stratify=y)
# X_train = np.concatenate((X_train, X_val))
# y_train = np.concatenate((y_train, y_val))

train_steps_per_epoch = X_train.shape[0] // batch_size
val_steps = X_val.shape[0] // batch_size

data_gen = ImageDataGenerator(rescale=1 / 255.)
training_generator = data_gen.flow(X_train, y_train, batch_size=batch_size)
validation_generator = data_gen.flow(X_val, y_val, batch_size=batch_size)

history = model.fit_generator(training_generator,
							  steps_per_epoch=train_steps_per_epoch,
							  epochs=epochs,
							  validation_data=validation_generator,
							  validation_steps=val_steps)

print(model.summary)
# from keras.utils import plot_model
# plot_model(conv_lstm_model, to_file='conv_lstm_model.png')


text:
#! /usr/bin/env python
#
# util_NRExtrudeOverlapGrid.py
#
#   This takes the whole NR catalog (aligned spin only!) and extrudes it in mtot (only)
#   Comparisons for cutoff are done ONLY using the 'approx' in LAL, for convenience
#
# EXAMPLES
#
#     # grid targeted to a single simulation
#        python util_NRExtrudeOverlapGrid.py --inj mdc.xml.gz  --group Sequence-SXS-All --param 1 --event 0 --verbose --skip-overlap
#     # extend the grid, to include the rest of the simlations
#        python util_NRExtrudeOverlapGrid.py --inj overlap-grid.xml.gz  --group Sequence-SXS-All  --event 0 --verbose --skip-overlap --insert-missing-spokes
#
# WARNINGS: i
#   xi_factor: meaning of grid is CHANGED from original interpretation
#

import argparse
import sys
import numpy as np
import RIFT.lalsimutils as lalsimutils
import lalsimulation as lalsim
import lalframe
import lal
import functools

import RIFT.physics.effectiveFisher  as eff   # for the mesh grid generation
import RIFT.physics.PrecessingFisherMatrix   as pcf   # Superior tools to perform overlaps. Will need to standardize with Evans' approach in effectiveFisher.py

from multiprocessing import Pool
try:
    import os
    n_threads = int(os.environ['OMP_NUM_THREADS'])
    print(" Pool size : ", n_threads)
except:
    n_threads=1
    print(" - No multiprocessing - ")

try:
	import NRWaveformCatalogManager3 as nrwf
	hasNR =True
except:
	hasNR=False
try:
    hasEOB=True
    import EOBTidalExternal as eobwf
except:
    hasEOB=False

###
### Load options
###

parser = argparse.ArgumentParser()
# Parameters
parser.add_argument("--mtot-range",default='[50,110]')
parser.add_argument("--mc-range",default=None, help='If specified, overrides the total mass range [24,38]')
parser.add_argument("--grid-cartesian-npts",default=30)
parser.add_argument("--group",default=None)
parser.add_argument("--param", action='append', help='Explicit list of parameters to use')
parser.add_argument("--insert-missing-spokes", action='store_true')
parser.add_argument("--aligned-only",action='store_true')
parser.add_argument("--eta-range",default='[0.1,0.25]')
parser.add_argument("--mass-xi-factor",default=0.0,type=float, help="The mass ranges are assumed to apply at ZERO SPIN. For other values of xi, the mass ranges map to m_{used} = m(1+xi *xi_factor+(1/4-eta)*eta_factor).  Note that to be stable, xi_factor<1. Default value 0.6, based on relevant mass region")
parser.add_argument("--mass-eta-factor",default=0,type=float, help="The mass ranges are assumed to apply at ZERO SPIN. For other values of xi, the mass ranges map to m_{used} = m(1+xi *xi_factor+(1/4-eta)*eta_factor).  Note that to be stable, xi_factor<1. Default value 0.6, based on relevant mass region")
# Cutoff options
parser.add_argument("--skip-overlap",action='store_true', help="If true, the grid is generated without actually performing overlaps. Very helpful if the grid is just in mtot, for the purposes of reproducing a specific NR simulation")
parser.add_argument("--match-value", type=float, default=0.01, help="Use this as the minimum match value. Default is 0.01 (i.e., keep almost everything)")
# Overlap options
parser.add_argument("--fisher-psd",type=str,default="SimNoisePSDaLIGOZeroDetHighPower",help="psd name ('eval'). lalsim., lalsim.SimNoisePSDaLIGOZeroDetHighPower, lalsimutils.Wrapper_AdvLIGOPsd, .SimNoisePSDiLIGOSRD... ")
parser.add_argument("--psd-file",  help="File name for PSD (assumed hanford). Overrides --fisher-psd if provided")
parser.add_argument("--srate",type=int,default=16384,help="Sampling rate")
parser.add_argument("--seglen", type=float,default=64., help="Default window size for processing. Short for NR waveforms")
parser.add_argument("--fref",type=float,default=0.);
# External grid
parser.add_argument("--external-grid-xml", default=None,help="Inspiral XML file (injection form) for alternate grid")
parser.add_argument("--external-grid-txt", default=None, help="Cartesian grid. Must provide parameter names in header. Exactly like output of code. Last column not used.")
# Base point
parser.add_argument("--inj", dest='inj', default=None,help="inspiral XML file containing the base point (OR containing the reference XML grid to check against a database)")
parser.add_argument("--event",type=int, dest="event_id", default=None,help="event ID of injection XML to use.")
parser.add_argument("--fmin", default=10,type=float,help="Mininmum frequency in Hz, default is 40Hz to make short enough waveforms. Focus will be iLIGO to keep comutations short")
parser.add_argument("--require-fmin-above-NR-start",action='store_true')
parser.add_argument("--fmax",default=2000,type=float,help="Maximum frequency in Hz, used for PSD integral.")
parser.add_argument("--mass1", default=35,type=float,help="Mass in solar masses")  # 150 turns out to be ok for Healy et al sims
parser.add_argument("--mass2", default=35,type=float,help="Mass in solar masses")
parser.add_argument("--s1z", default=0.1,type=float,help="Spin1z")
#parser.add_argument("--lambda1",default=590,type=float)
#parser.add_argument("--lambda2", default=590,type=float)
parser.add_argument("--eff-lambda", type=float, help="Value of effective tidal parameter. Optional, ignored if not given")
parser.add_argument("--deff-lambda", type=float, help="Value of second effective tidal parameter. Optional, ignored if not given")
parser.add_argument("--lmax", default=2, type=int)
parser.add_argument("--approx",type=str,default="SEOBNRv2")
# Output options
parser.add_argument("--fname", default="overlap-grid", help="Base output file for ascii text (.dat) and xml (.xml.gz)")
parser.add_argument("--verbose", action="store_true",default=False, help="Required to build post-frame-generating sanity-test plots")
parser.add_argument("--save-plots",default=False,action='store_true', help="Write plots to file (only useful for OSX, where interactive is default")
opts=  parser.parse_args()

xi_factor = opts.mass_xi_factor
eta_factor0 = opts.mass_eta_factor

if opts.verbose:
    True
    #lalsimutils.rosDebugMessagesContainer[0]=True   # enable error logging inside lalsimutils

###
### Define grid overlap functions
###   - Python's 'multiprocessing' module seems to cause process lock
###

Lmax = 2

def eval_overlap(grid,P_list, IP,indx):
#    if opts.verbose:
#        print " Evaluating for ", indx
    global Lmax
    global opts
    P2 = P_list[indx]
#    P2.print_params()
    line_out = []
    line_out = list(grid[indx])
    if not opts.skip_overlap:
        T_here = 1./IP.deltaF
        P2.deltaF=1./T_here
        hf2 = lalsimutils.complex_hoff(P2)
        nm2 = IP.norm(hf2);  hf2.data.data *= 1./nm2
        ip_val = IP.ip(hfBase,hf2)
        line_out.append(ip_val)
    else:
        line_out.append(-1)
    if opts.verbose:
            print(" Answer ", indx, line_out)
    return line_out

def evaluate_overlap_on_grid(hfbase,param_names, grid):
    # Validate grid is working: Create a loop and print for each one.
    # WARNING: Assumes grid for mass-unit variables hass mass units (!)
    global xi_factor
    P_list = []
    for line in grid:
        Pgrid = P.manual_copy()
        # Set attributes that are being changed as necessary, leaving all others fixed
        for indx in np.arange(len(param_names)):
                Pgrid.assign_param(param_names[indx], line[indx])   # DANGER: need to keep overlap values
        # Rescale mass parameters using the xi factor
        xi_now = Pgrid.extract_param('xi')
        eta_now = Pgrid.extract_param('eta')
            # THIS CAN BE NEGATIVE: PATHOLOGICAL
        Pgrid.m1 *= (1.+xi_now*xi_factor - eta_factor0*(0.25-eta_now))
        Pgrid.m2 *= (1.+xi_now*xi_factor - eta_factor0*(0.25-eta_now))
        P_list.append(Pgrid)
#    print "Length check", len(P_list), len(grid)
    ###
    ### Loop over grid and make overlaps : see effective fisher code for wrappers
    ###
    #  FIXME: More robust multiprocessing implementation -- very heavy!
#    p=Pool(n_threads)
    # PROBLEM: Pool code doesn't work in new configuration.
    grid_out = np.array(map(functools.partial(eval_overlap, grid, P_list,IP), np.arange(len(grid))))
    # Remove mass units at end
    for p in ['mc', 'm1', 'm2', 'mtot']:
        if p in param_names:
            indx = param_names.index(p)
            grid_out[:,indx] /= lal.MSUN_SI
    # Truncate grid so overlap with the base point is > opts.min_match. Make sure to CONSISTENTLY truncate all lists (e.g., the P_list)
    grid_out_new = []
    P_list_out_new = []
    for indx in np.arange(len(grid_out)):
      if P_list[indx].m1 >0 and P_list[indx].m2>0:   # reject insane cases where the masses were scaled to be <0
        if (opts.skip_overlap) or (grid_out[indx,-1] > opts.match_value):
            grid_out_new.append(grid_out[indx])
            P_list_out_new.append(P_list[indx])
    grid_out = np.array(grid_out_new)
    return grid_out, P_list_out_new

###
### Define base point
###

# Handle PSD
# FIXME: Change to getattr call, instead of 'eval'
eff_fisher_psd = lalsim.SimNoisePSDiLIGOSRD
if not opts.psd_file:
    #eff_fisher_psd = eval(opts.fisher_psd)
    eff_fisher_psd = getattr(lalsim, opts.fisher_psd)   # --fisher-psd SimNoisePSDaLIGOZeroDetHighPower   now
    analyticPSD_Q=True
else:
    eff_fisher_psd = lalsimutils.load_resample_and_clean_psd(opts.psd_file, 'H1', 1./opts.seglen)
    analyticPSD_Q = False

P=lalsimutils.ChooseWaveformParams()
if opts.inj:
    from glue.ligolw import lsctables, table, utils # check all are needed
    filename = opts.inj
    event = opts.event_id
    xmldoc = utils.load_filename(filename, verbose = True,contenthandler =lalsimutils.cthdler)
    sim_inspiral_table = table.get_table(xmldoc, lsctables.SimInspiralTable.tableName)
    P.copy_sim_inspiral(sim_inspiral_table[int(event)])
else:
    P.m1 = opts.mass1 *lal.MSUN_SI
    P.m2 = opts.mass2 *lal.MSUN_SI
    P.s1z = opts.s1z
    P.dist = 150*1e6*lal.PC_SI
    if opts.eff_lambda and Psig:
        lambda1, lambda2 = 0, 0
        if opts.eff_lambda is not None:
            lambda1, lambda2 = lalsimutils.tidal_lambda_from_tilde(m1, m2, opts.eff_lambda, opts.deff_lambda or 0)
            Psig.lambda1 = lambda1
            Psig.lambda2 = lambda2

    P.fmin=opts.fmin   # Just for comparison!  Obviously only good for iLIGO
    P.ampO=-1  # include 'full physics'
    if opts.approx:
        P.approx = lalsim.GetApproximantFromString(opts.approx)
        if not (P.approx in [lalsim.TaylorT1,lalsim.TaylorT2, lalsim.TaylorT3, lalsim.TaylorT4]):
            # Do not use tidal parameters in approximant which does not implement them
            print(" Do not use tidal parameters in approximant which does not implement them ")
            P.lambda1 = 0
            P.lambda2 = 0
    else:
        P.approx = lalsim.GetApproximantFromString("TaylorT4")
P.deltaT=1./16384
P.taper = lalsim.SIM_INSPIRAL_TAPER_START
P.deltaF = 1./opts.seglen #lalsimutils.findDeltaF(P)
P.fref = opts.fref
P.print_params()

print("    -------INTERFACE ------")
hfBase =None
IP=None
if not opts.skip_overlap:
    hfBase = lalsimutils.complex_hoff(P)
    IP = lalsimutils.CreateCompatibleComplexOverlap(hfBase,analyticPSD_Q=analyticPSD_Q,psd=eff_fisher_psd,fMax=opts.fmax)
    nmBase = IP.norm(hfBase)
    hfBase.data.data *= 1./nmBase
    if opts.verbose:
        print(" ------  SIGNAL DURATION ----- ")
        print(hfBase.data.length*P.deltaT)

###
### If we are testing an existing grid to see if simulations are missing
###

# Step 1: Load in the grid, identify NR simulations

P_list_nr = []
gp_list = []
if opts.inj and opts.insert_missing_spokes:
    import spokes
    sdHere = spokes.LoadSpokeXML(opts.inj)
    for key in sdHere:
        P_here =sdHere[key][0]
        P_list_nr.append(P_here)
        compare_dict = {}
        compare_dict['q'] = P_here.m2/P_here.m1 # Need to match the template parameter. NOTE: VERY IMPORTANT that P is updated with the event params
        compare_dict['s1z'] = P_here.s1z
        compare_dict['s1x'] = P_here.s1x
        compare_dict['s1y'] = P_here.s1y
        compare_dict['s2z'] = P_here.s2z
        compare_dict['s2x'] = P_here.s2x
        compare_dict['s2y'] = P_here.s2y
        good_sim_list = nrwf.NRSimulationLookup(compare_dict,valid_groups=[opts.group])
        if len(good_sim_list)< 1:
                        print(" ------- NO MATCHING SIMULATIONS FOUND ----- ")
                        import sys
                        sys.exit(0)
                        print(" Identified set of matching NR simulations ", good_sim_list)
        try:
                        print("   Attempting to pick longest simulation matching  the simulation  ")
                        MOmega0  = 1
                        good_sim = None
                        for key in good_sim_list:
                                print(key, nrwf.internal_EstimatePeakL2M2Emission[key[0]][key[1]])
                                if nrwf.internal_WaveformMetadata[key[0]][key[1]]['Momega0'] < MOmega0:
                                        good_sim = key
                                        MOmega0 = nrwf.internal_WaveformMetadata[key[0]][key[1]]['Momega0']
                                print(" Picked  ",key,  " with MOmega0 ", MOmega0, " and peak duration ", nrwf.internal_EstimatePeakL2M2Emission[key[0]][key[1]])
        except:
            good_sim  = good_sim_list[0] # pick the first one.  Note we will want to reduce /downselect the lookup process
            group = good_sim[0]
            param = good_sim[1]
            gp_list.append( (group,param))

###
### Load in the NR simulation array metadata
###
P_list_NR = []
omega_list_NR=[]

glist = []
if opts.group:
    glist = [opts.group]
else:
    glist = nrwf.internal_ParametersAvailable.keys()

eta_range = eval(opts.eta_range)

print(opts.group, opts.param)
for group in glist:
  if not opts.param:
    for param in nrwf.internal_ParametersAvailable[group]:
     print("  -> ", group, param)
     try:
        wfP = nrwf.WaveformModeCatalog(group,param,metadata_only=True)
        wfP.P.deltaT = P.deltaT
        wfP.P.deltaF = P.deltaF
        wfP.P.fmin = P.fmin
#        wfP.P.print_params()
        # Add parameters. Because we will compare with SEOB, we need an ALIGNED waveform, so we fake it
        if wfP.P.extract_param('eta') >= eta_range[0] and wfP.P.extract_param('eta')<=eta_range[1] and (not np.isnan(wfP.P.s1z) and not np.isnan(wfP.P.s2z)):
            if (not opts.skip_overlap) and wfP.P.SoftAlignedQ():
                print(" Adding aligned sim ", group, param,  " and changing spin values, so we can compute overlaps ")
                wfP.P.approx = lalsim.GetApproximantFromString(opts.approx)  # Make approx consistent and sane
                wfP.P.m2 *= 1. - 1e-6  # Prevent failure for exactly equal!
            # Satisfy error checking condition for lal
                wfP.P.s1x = 0
                wfP.P.s2x = 0
                wfP.P.s1y = 0
                wfP.P.s2y = 0
                P_list_NR = P_list_NR + [wfP.P]
                try:
                    omega_list_NR += [nrwf.internal_WaveformMetadata[group][param]["Momega0"]]
                except:
                    omega_list_NR+=[-1]
            elif opts.skip_overlap:
                if not opts.aligned_only:
                    print(" Adding generic sim; for layout only ", group, param)
                    P_list_NR = P_list_NR + [wfP.P]
                    try:
                        omega_list_NR += [nrwf.internal_WaveformMetadata[group][param]["Momega0"]]
                    except:
                        omega_list_NR+=[-1]
                elif opts.aligned_only and  wfP.P.SoftAlignedQ():
                    print(" Adding aligned spin simulation; for layout only", group, param, " and fixing transverse spins accordingly")
                    wfP.P.s1x = 0
                    wfP.P.s2x = 0
                    wfP.P.s1y = 0
                    wfP.P.s2y = 0
                    P_list_NR = P_list_NR + [wfP.P]
                    try:
                        omega_list_NR += [nrwf.internal_WaveformMetadata[group][param]["Momega0"]]
                    except:
                        omega_list_NR+=[-1]

            else:
                print(" Skipping non-aligned simulation because overlaps active (=SEOBNRv2 comparison usually)", group, param)
#                wfP.P.print_params()
#                print nrwf.internal_WaveformMetadata[group][param]
     except:
	print(" Failed to add ", group, param)

  else: # target case if a single group and parameter sequence are specified
        print("Looping over list ", opts.param)
        for paramKey in opts.param:
            print(" Adding specific simulation ", opts.group, paramKey)
            if nrwf.internal_ParametersAreExpressions[group]:
                param = eval(paramKey)
            else:
                param = paramKey
        wfP = nrwf.WaveformModeCatalog(group,param,metadata_only=True)
        wfP.P.deltaT = P.deltaT
        wfP.P.deltaF = P.deltaF
        wfP.P.fmin = P.fmin
        P_list_NR = P_list_NR + [wfP.P]
        try:
            omega_list_NR += [nrwf.internal_WaveformMetadata[group][param]["Momega0"]]
        except:
            omega_list_NR+=[-1]

if len(P_list_NR)<1:
    print(" No simulations")
    sys.exit(0)

###
### Define parameter ranges to be changed
###

template_min_freq = opts.fmin
ip_min_freq = opts.fmin

###
### Lay out grid, currently CARTESIAN.   OPTIONAL: Load grid from file
###

# For now, we just extrude in these parameters
if not opts.skip_overlap:  # aligned only!  Compare to SEOB
    if opts.mc_range:
        param_names = ['mc', 'eta','s1z','s2z']
    else:
        param_names = ['mtot', 'q','s1z','s2z']
else:
    if opts.mc_range:
        param_names = ['mc', 'eta', 's1x', 's1y', 's1z', 's2x', 's2y', 's2z']
    else:
        param_names = ['mtot', 'q', 's1x', 's1y', 's1z', 's2x', 's2y', 's2z']

mass_range =[]
if opts.mc_range:
    mass_range = np.array(eval(opts.mc_range))*lal.MSUN_SI
else:
    mass_range = np.array(eval(opts.mtot_range))*lal.MSUN_SI
mass_grid =np.linspace( mass_range[0],mass_range[1],opts.grid_cartesian_npts)

# Loop over simulations and mass grid
grid = []
indx=-1
print("Length check", len(P_list_NR), len(omega_list_NR))
for P in P_list_NR:
    indx+=1
    Momega0_here = omega_list_NR[indx]
#    P.print_params()
    eta0 = P.extract_param('eta')   # prevent drift in mass ratio across the set if I change mc
    for M in mass_grid:
        if opts.mc_range:
            P.assign_param('mc', M)
            P.assign_param('eta', eta0)  # prevent drift in mass ratio
        else:
            P.assign_param('mtot', M)
        newline = []
        for param in param_names:
            newline = newline + [P.extract_param(param)]
#        print newline
        fstart_NR = Momega0_here/((P.extract_param('mtot')/(lal.MSUN_SI) ) *lalsimutils.MsunInSec*(np.pi))
        if   fstart_NR < P.fmin or not opts.require_fmin_above_NR_start :
            # only add mass point if a valid mass choice, given NR starting frequency.
            grid = grid+[ newline]

print(" ---- DONE WITH GRID SETUP --- ")
print(" grid points # = " ,len(grid))

grid_out, P_list = evaluate_overlap_on_grid(hfBase, param_names, grid)
if len(grid_out)==0:
    print(" No points survive....")

###
### Optional: Write grid to XML file (ONLY if using cutoff option)
###
lalsimutils.ChooseWaveformParams_array_to_xml(P_list, fname=opts.fname, fref=P.fref)

###
### Write output to text file:  p1 p2 p3 ... overlap, only including named params
###
headline = ' '.join(param_names + ['ip'])
if int(np.version.short_version.split('.')[1])>=7:  # http://docs.scipy.org/doc/numpy/reference/generated/numpy.savetxt.html
    np.savetxt(opts.fname+".dat", grid_out, header=headline)
else:
    np.savetxt(opts.fname+".dat", grid_out)   #

###
### Optional: Scatterplot
###
if opts.verbose and len(param_names)==1 and len(grid_out)>0:
    import matplotlib.pyplot as plt
    fig = plt.figure()
    plt.plot(grid_out[:,0], grid_out[:,1])
    plt.show()

if opts.verbose and len(param_names)==2:
    import matplotlib.pyplot as plt
    from mpl_toolkits.mplot3d import Axes3D
    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')
    ax.scatter(grid_out[:,0], grid_out[:,1], grid_out[:,2])
    plt.show()

print(" ---- DONE ----")

max_stars_repo_name: spfanning/research-projects-RIT

max_issues_repo_name: spfanning/research-projects-RIT

max_forks_repo_name: spfanning/research-projects-RIT


text:
import pickle
import operator
from collections import Counter

from FormulaRetrieval import getFormulae

from arXivDocs2Vec import docs2vec
from arXivDocs2tfidf import docs2tfidf

from fuzzywuzzy import fuzz

import numpy as np
from sklearn.neighbors import NearestNeighbors

# set file paths
#basePath = ""
#basePath = "F:\\SigMathLing\\"
basePath = "F:\\arXiv\\"
inputPath = basePath + "formulae\\mathmlben\\"
outputPath = basePath + "formulae\\mathmlben\\"

# load mmlben arXiv matches results
# with open(outputPath + "results_mmlben_arXiv.pkl", 'rb') as f:
#     matches = pickle.load(f)
with open(outputPath + "matches_EFE_arXiv.pkl", 'rb') as f:
    matches = pickle.load(f)

for formula in matches.items():
    print(formula[0])

# load formula catalog mmlben
# with open(outputPath + "formula_catalog_mmlben.pkl",'rb') as f:
#     formula_catalog_mmlben = pickle.load(f)
# load formula catalog arXiv all
with open(outputPath + "formula_catalog_all.pkl",'rb') as f:
    formula_catalog_all = pickle.load(f)

# formula catalog arXiv
# create and save
# with open(outputPath + "formula_catalog_backup.pkl",'rb') as f:
#     formula_catalog_backup = pickle.load(f)
# del formula_catalog_backup['operator_catalog']
# del formula_catalog_backup['identifier_catalog']
# formula_catalog_arXiv = {}
# for formula in formula_catalog_backup.items():
#     try:
#         print(formula[1]['TeX'])
#         formula_catalog_arXiv[formula[1]['TeX']] = {}
#         formula_catalog_arXiv[formula[1]['TeX']]['operators'] = []
#         formula_catalog_arXiv[formula[1]['TeX']]['identifiers'] = []
#         for operator in formula[1]['operators'].values():
#             formula_catalog_arXiv[formula[1]['TeX']]['operators'].append(operator)
#         for identifier in formula[1]['identifiers'].values():
#             formula_catalog_arXiv[formula[1]['TeX']]['identifiers'].append(identifier)
#     except:
#         pass
# with open(outputPath + "formula_catalog_arXiv.pkl",'wb') as f:
#     pickle.dump(formula_catalog_arXiv,f)
# load
# with open(outputPath + "formula_catalog_arXiv.pkl",'rb') as f:
#     formula_catalog_arXiv = pickle.load(f)

# formula catalog mmlben
# create and save
# with open(outputPath + "mathmlben.html", 'r', encoding="utf8") as f:
#     formula_catalog_mmlben = getFormulae(doc_str=f.read(),mode='math')
# with open(outputPath + "formula_catalog_mmlben.pkl",'wb') as f:
# #     pickle.dump(formula_catalog_mmlben,f)
# load
# with open(outputPath + "formula_catalog_mmlben.pkl",'rb') as f:
#     formula_catalog_mmlben = pickle.load(f)

# create formula labs and data list
#formulaLabs = []
#formulaData = []
# bring formula catalog operators and identifiers in data array
# def catalog_to_data(formula_catalog):
#     for formula in formula_catalog.items():
#         # data needs to be a string per formula
#         Data = ""
#         for operator in formula[1]['operators']:
#             Data += operator + " "
#         for identifier in formula[1]['identifiers']:
#             Data += identifier + " "
#         # cut off last whitespace at the end
#         formulaData.append(Data[:-1])
#         # tex string as label
#         formulaLabs.append(formula[0])
# insert formula catalog mmlben
#catalog_to_data(formula_catalog_mmlben)
# insert formula catalog arXiv
#catalog_to_data(formula_catalog_arXiv)

# save data and labels
# with open(outputPath + "formulaData.pkl","wb") as f:
#     pickle.dump(formulaData,f)
# with open(outputPath + "formulaLabs.pkl","wb") as f:
#     pickle.dump(formulaLabs,f)

# build and save doc2Vec vectors
# model,formulaVecs_doc2vec = docs2vec(formulaData,formulaLabs)
# with open(outputPath + "formulaVecs_doc2vec.pkl","wb") as f:
#     pickle.dump(formulaVecs_doc2vec,f)
# # build and save tf-idf vectors
# formulaVecs_tfidf = docs2tfidf(formulaData)
# with open(outputPath + "formulaVecs_tfidf.pkl","wb") as f:
#     pickle.dump(formulaVecs_tfidf,f)

# load data and labels
# with open(outputPath + "formulaData.pkl","rb") as f:
#     formulaData = pickle.load(f)
# with open(outputPath + "formulaLabs.pkl","rb") as f:
#     formulaLabs = pickle.load(f)
# load doc2vec and tf-idf vectors
# with open(outputPath + "formulaVecs_doc2vec.pkl","rb") as f:
#     formulaVecs_doc2vec = pickle.load(f)
# with open(outputPath + "formulaVecs_tfidf.pkl","rb") as f:
#     formulaVecs_tfidf = pickle.load(f)

# load FCR candidates
# with open(outputPath + "results_fuzzy.pkl","rb") as f:
#     results_fuzzy = pickle.load(f)
# with open(outputPath + "results_seq.pkl","rb") as f:
#     results_seq = pickle.load(f)
# with open(outputPath + "results_set.pkl","rb") as f:
#     results_set = pickle.load(f)
# with open(outputPath + "results_50.pkl","rb") as f:
#     results_50 = pickle.load(f)

# mmlben formulaLab nr.s neq GoldIDs
# selection = [24,219,275,285,294,295,300,301,302,303]
# formula_catalog_selection = {}
# for nr in selection:
#     formula_catalog_selection[formulaLabs[nr]] = {}
#     formula_catalog_selection[formulaLabs[nr]]['operators'] = set(formula_catalog_mmlben[formulaLabs[nr]]['operators'])
#     formula_catalog_selection[formulaLabs[nr]]['identifiers'] = set(formula_catalog_mmlben[formulaLabs[nr]]['identifiers'])
#selection = list(range(0,302))

# without encoding

# tex
# fuzzy formula tex string search
# results_fuzzy = {}
# for nr in selection:
#     formulaQuery = formulaLabs[nr]
#     #formulaQuery = "'\\frac {1}{c^2} \\frac{\\partial^2}{\\partial t^2} \\psi - \\nabla^2 \\psi + \\frac {m^2 c^2}{\\hbar^2} \\psi = 0'"
#     results = {}
#     i = 0
#     for formulaLab in formulaLabs:
#         print(str(round(np.multiply(np.divide(i,len(formulaLabs)),100),1)) + " %")
#         try:
#             fpr = fuzz.partial_ratio(formulaQuery,formulaLab)
#             if fpr > 80 and len(formulaLabs) > 3:
#                 results[formulaLab] = fpr
#         except:
#             pass
#         i += 1
#     results_fuzzy[formulaLabs[nr]] = sorted(results.items(),key=operator.itemgetter(1),reverse=True)

# mml
# fuzzy formula parts seq string search
# results_seq = {}
# for nr in selection:
#     formulaQuery = formulaLabs[nr]
#     #formulaQuery = "'\\frac {1}{c^2} \\frac{\\partial^2}{\\partial t^2} \\psi - \\nabla^2 \\psi + \\frac {m^2 c^2}{\\hbar^2} \\psi = 0'"
#     results = {}
#     i = 0
#     for formulaLab in formulaLabs:
#         print(str(round(np.multiply(np.divide(i,len(formulaLabs)),100),1)) + " %")
#         try:
#             fpr = fuzz.partial_ratio(formulaQuery,formulaData[i])
#             if fpr > 80 and len(formulaLabs) > 3:
#                 results[formulaLab] = fpr
#         except:
#             pass
#         i += 1
#     results_seq[formulaLabs[nr]] = sorted(results.items(),key=operator.itemgetter(1),reverse=True)

# search for mathmlben formulae by parts (op,id) set
#results_set = {}
# #for nr in selection:
# for nr in range(0,50):
#     print(nr)
#     formulaQuery = set(formulaData[nr].split())
#     #formulaQuery = set(['∂','ψ','t','V','∇','+','>','m','c','ℏ'])
#     #formulaQuery = set(['G','μ','ν','T','κ','π','c'])
#     results = {}
#     for i in range(0,len(formulaData)):
#         results[formulaLabs[i]] = len(formulaQuery.intersection(set(formulaData[i].split())))
#     results_set[formulaLabs[nr]] = sorted(results.items(), key=operator.itemgetter(1),reverse=True)

# FCR here

# results_set = {}
# for formula_query in formula_catalog_selection.items():
#     print("Query: " + formula_query[0])
#     results = {}
#     for formula_index in formula_catalog_all.items():
#         try:
#             overlap = len(formula_query[1]['identifiers'].intersection(formula_index[1]['identifiers']))# + len(formula_query[1]['operators'].intersection(formula_index[1]['operators']))
#             if overlap > 0:
#                 results[formula_index[0]] = overlap
#                 print(formula_index[0] + ": " + str(overlap))
#         except:
#             pass
#     results_set[formula_query[0]] = sorted(results.items(), key=operator.itemgetter(1),reverse=True)

# with encoding
# set encoding
#encoding = formulaVecs_doc2vec
#encoding = formulaVecs_tfidf.toarray()

# cosine similarity neighbors
# def cos_mat_mult (Vec, Mat):
#     Sim = []
#     for Vec2 in Mat:
#         Sim.append(float(np.inner(Vec, Vec2)))
#     return Sim

# euclidean distance neighbors
# def eucl_dist_vec (Vec, Mat):
#     Sim = []
#     for Vec2 in Mat:
#         Sim.append(float(np.linalg.norm(Vec-Vec2)))
#     return Sim

# search for mmlben formulae
# search_id = 1
# search_results = {}
# calculate cosine similarity with all vectors
#Sim = cos_mat_mult(encoding[idx],encoding)
# calculate euclidean distance to all vectors
#Sim = eucl_dist_vec(encoding[search_id],encoding)
# sort and label results
# for i in range(0,len(Sim)):
#     search_results[formulaLabs[i]] = Sim[i]
# search_results = sorted(search_results.items(),key=operator.itemgetter(1),reverse=True)

# k nearest neighbors
# nbrs = NearestNeighbors(n_neighbors=3, algorithm='auto').fit(encoding)
# distances, indices = nbrs.kneighbors(encoding)
#
# with open(outputPath + "kNN_indices.pkl",'wb') as f:
#     pickle.dump(indices,f)

#Query equation by matching identifers

# Klein-Gordon
#formulaQuery = set(['∂','ψ','t','V','∇','+','>','m','c','ℏ'])
# Schrödinger
#formulaQuery = set(['ψ','t','m','c','ℏ'])
# Einstein
formulaQuery = set(['G', 'R', 'T', 'g', 'μ', 'ν', 'Λ', 'c', 'κ'])

matches = {}
arXiv_counter = 0
for formula_index in formula_catalog_all.items():
    arXiv_counter += 1
    print(arXiv_counter)
    # nr_identifiers = len(formula_index[1]['identifiers'])
    # same_identifiers = len(formula_index[1]['identifiers'].difference(formulaQuery)) == 0
    overlap = len(formulaQuery.intersection(formula_index[1]['identifiers']))
    if overlap > 5:
        print(formula_index[0] + ": " + str(formula_index[1]))
        matches[formula_index[0]] = str(formula_index[1])

with open(outputPath + "matches_EFE_arXiv.pkl", 'wb') as f:
    pickle.dump(matches, f)

# Query mathmlben formulae
# matches = {}
# mmlben_counter = 0
# arXiv_counter = 0
# for formula_query in formula_catalog_mmlben.items():
#     mmlben_counter += 1
#     print(formula_query[0] + ": " + str(formula_query[1]))
#     for formula_index in formula_catalog_all.items():
#         arXiv_counter += 1
#         print(str(mmlben_counter) + ": " + arXiv_counter)
#         nr_operators = len(formula_index[1]['operators'])
#         nr_identifiers = len(formula_index[1]['identifiers'])
#         same_operators = len(formula_index[1]['operators'].difference(formula_query[1]['operators'])) == 0
#         same_identifiers = len(formula_index[1]['identifiers'].difference(formula_query[1]['identifiers'])) == 0
#         if same_operators and same_identifiers and nr_operators > 0 and nr_identifiers > 1:
#             print(formula_index[0] + ": " + str(formula_index[1]))
#             matches[formula_index[0]] = str(formula_index[1])

# with open(outputPath + "results_id_mmlben_arXiv.pkl",'wb') as f:
#      pickle.dump(results_set,f)

print("end")

max_stars_repo_name: philsMINT/formula-concept-retrieval

max_issues_repo_name: philsMINT/formula-concept-retrieval

max_forks_repo_name: philsMINT/formula-concept-retrieval


text:
#!/usr/bin/env python3

import argparse

import pandas as pd
import numpy as np
from sklearn.neighbors import KDTree

def compute_neighbors(tsne_file):
    # Parse t-SNE coordinates file
    df1 = pd.read_csv(tsne_file, sep='\t')

    # Build a KDTree with the coordinates
    tree = KDTree(df1[['t1', 't2']].values)

    # Pull user input samples to query
    inputs = df1[df1['projects'] == 'input']
    inputs = inputs.reset_index()
    inputs_pos = inputs[['t1', 't2']].values

    summary = "{0}\t{1}\t{2}\t{3}\t{4}\t{5}\t{6}\t{7}\t{8}\n".format("sample", "top disease code from n5", "support for top disease code from n5", 'n5 neighbors', 'n5 disease counts', "top disease code from n10", "support for top disease code from n10", 'n10 neighbors', 'n10 disease counts')
    # For each input sample, compute 5 nearest neighbors
    for index, row in inputs.iterrows():
        query = np.reshape(row[['t1','t2']].values, (1,2))

        dist, ind = tree.query(query, k=6)

        disease = {}
        neighbors = ""

        for val in ind[0,:]:
            if not df1.loc[val, 'samples'] == row['samples']:
                if df1.loc[val, 'classes'] in disease:
                    disease[df1.loc[val, 'classes']] += 1
                else:
                    disease[df1.loc[val, 'classes']] = 1
                if neighbors == "":
                    neighbors += df1.loc[val, 'samples'] + "," + df1.loc[val, 'classes']
                else:
                    neighbors += ";" + df1.loc[val, 'samples'] + "," + df1.loc[val, 'classes']

        counts=""
        top=""
        top_count=""
        for key, value in sorted(disease.items(), key=lambda item: item[1], reverse=True):
            if counts == "":
                top=key
                top_count=value
                counts += "{0}:{1}".format(key, value)
            else:
                counts += ";{0}:{1}".format(key, value)

        dist_n10, ind_n10 = tree.query(query, k=11)
        disease_n10 = {}
        neighbors_n10 = ""

        for val in ind_n10[0,:]:
            if not df1.loc[val, 'samples'] == row['samples']:
                if df1.loc[val, 'classes'] in disease_n10:
                    disease_n10[df1.loc[val, 'classes']] += 1
                else:
                    disease_n10[df1.loc[val, 'classes']] = 1
                if neighbors_n10 == "":
                    neighbors_n10 += df1.loc[val, 'samples'] + "," + df1.loc[val, 'classes']
                else:
                    neighbors_n10 += ";" + df1.loc[val, 'samples'] + "," + df1.loc[val, 'classes']

        counts_n10=""
        top_n10=""
        top_count_n10=""
        for key, value in sorted(disease_n10.items(), key=lambda item: item[1], reverse=True):
            if counts_n10 == "":
                top_n10=key
                top_count_n10=value
                counts_n10 += "{0}:{1}".format(key, value)
            else:
                counts_n10 += ";{0}:{1}".format(key, value)

        summary += "{0}\t{1}\t{2}\t{3}\t{4}\t{5}\t{6}\t{7}\t{8}\n".format(row['samples'], top, top_count, neighbors, counts, top_n10, top_count_n10, neighbors_n10, counts_n10)

    return summary

def main():
    parser = argparse.ArgumentParser()
    parser.add_argument("tsne", help="t-SNE coordinates file to process")
    parser.add_argument("output_file", default="neighbors.tsv", help="Output file to write neighborhood information")
    args = parser.parse_args()

    neighbors = compute_neighbors(args.tsne)

    with open(args.output_file, "w") as file:
        file.write(neighbors)

if __name__ == "__main__":
    main()

max_stars_repo_name: stjudecloud/expression-classification

max_issues_repo_name: stjudecloud/expression-classification

max_forks_repo_name: stjudecloud/expression-classification


text:
import numpy as np
from cost_functions import trajectory_cost_fn
import time

class Controller():
    def __init__(self):
        pass

    # Get the appropriate action(s) for this state(s)
    def get_action(self, state):
        pass

class RandomController(Controller):
    def __init__(self, env):
        """ YOUR CODE HERE """
        self.env = env
        self.ac_dim = env.action_space.shape[0]

    def get_action(self, state):
        """ YOUR CODE HERE """
        """ Your code should randomly sample an action uniformly from the action space """
        return np.random.uniform(low=-1.0, high=1.0, size=self.ac_dim)

class MPCcontroller(Controller):
    """ Controller built using the MPC method outlined in https://arxiv.org/abs/1708.02596 """
    def __init__(self,
                 env,
                 dyn_model,
                 horizon=5,
                 cost_fn=None,
                 num_simulated_paths=10,
                 ):
        self.env = env
        self.dyn_model = dyn_model
        self.horizon = horizon
        self.cost_fn = cost_fn
        self.num_simulated_paths = num_simulated_paths
        self.ac_dim = env.action_space.shape[0]

    def get_action(self, state):
        """ YOUR CODE HERE """
        """ Note: be careful to batch your simulations through the model for speed """
        obs = np.resize(state, (self.num_simulated_paths, len(state)))

        observations, next_observations = [], []
        actions = np.random.uniform(low=-1.0, high=1.0, size=(self.horizon, self.num_simulated_paths, self.ac_dim))

        for i in range(self.horizon):
            observations.append(obs)
            obs = self.dyn_model.predict(obs, actions[i,::])
            next_observations.append(obs)

        observations = np.array(observations)
        next_observations = np.array(next_observations)
        cost = trajectory_cost_fn(self.cost_fn, observations, actions, next_observations)
        min_idx = np.argmin(cost)
        return actions[0,min_idx,:]


text:
import os
import pandas as pd
import numpy as np
from ixtlilton_tools._private_tools.exceptions import UserDoesNotExist, NoAdminRights
import subprocess
import pathlib

def remove(username):

    from ixtlilton_tools.admin import running_with_admin_rights
    from ixtlilton_tools.admin.user import exists as user_exists
    from ixtlilton_tools.admin.directory.user import scratch, work, home

    if not running_with_admin_rights():
        raise NoAdminRights()

    if not user_exists(username):
        raise UserDoesNotExist(username=username)

    ## removing user

    subprocess.run(['userdel', username])

    ## removing mailbox file

    path=pathlib.Path('/var/spool/mail/'+username)
    if path.exists():
        os.remove(path)

    ## removing directories if empty

    scratch.remove(username)
    work.remove(username)
    home.local.opt.modulefiles.remove(username)
    home.local.opt.apps.remove(username)
    home.local.opt.remove(username)
    home.local.bin.remove(username)
    home.local.etc.remove(username)
    home.local.lib.remove(username)
    home.local.lib32.remove(username)
    home.local.lib64.remove(username)
    home.local.include.remove(username)
    home.local.src.remove(username)
    home.local.share.remove(username)
    home.local.man.remove(username)
    home.local.remove(username)
    home.remove(username)

max_stars_repo_name: uibcdf/Ixtlilton

max_issues_repo_name: uibcdf/Ixtlilton

max_forks_repo_name: uibcdf/Ixtlilton


text:
# -*- coding: utf-8 -*-

"""
ソルバーのクラス
"""

import numpy as np
import copy
import time
from src.const import Const

class Solver:

  ## public function ##

  def __init__(self):
    # self.board = board
    self.select_result_log = []

  # メイン関数
  def solve(self, board):
    while True:
      # ロジカルに進められるルーティン
      update_flg = self.update_routine_logical(board)

      if board.is_complete():
        print("===== COMPLETE! =====")
        break

      if not(update_flg):
        # 現状のデータをダンプする
        board.store_board_data()
        # 探索するセルを選択する
        select_result = self.select_search_cell(board)
        select_pos, select_val = select_result
        # 探索履歴を保持する
        self.select_result_log.append(select_result)
        board.update_cell(select_pos, select_val)
        board.update_candidate_list(select_pos, select_val)

      if board.is_wrong():
        board.restore_board_data()
        select_result = self.select_result_log.pop()
        select_pos, select_val = select_result
        # 候補リストから探索した結果を除外する
        board.remove_from_candidate_list(select_pos, select_val)

  ## private function ##

  # 論理的な手法で更新するルーチン
  def update_routine_logical(self, board):
    update_flg_1 = self.update_data_with_candidate_list_cell(board)
    update_flg_2 = self.update_data_with_candidate_list_area(board)

    update_flg_3 = board.update_all_candidate_list()

    if not(update_flg_1 or update_flg_2 or update_flg_3):
      update_flg_4 = board.update_candidate_list_2()
      if not(update_flg_4):
        return False

    return True

  # 候補リストのうち、1つになった値をセルに入れる
  def update_data_with_candidate_list_cell(self, board):
    # 更新があったかどうか
    update_flg = False
    for key, value in board.candidate_list[Const.CELL_KEY_NAME].items():
      if len(value) == 1:
        update_flg = True
        val = value.pop()
        board.update_cell(key, val)
    return update_flg

  # 特定の縦・横・矩形において、候補リストから各値の出現回数を取得する
  def update_data_with_candidate_list_area(self, board):
    update_flg = False
    # 領域(行・列・矩形)内における各値が入る可能性がある候補セルのリストがただ1つであればそのセルに値を入れる
    for mode, value in board.candidate_list[Const.AREA_KEY_NAME].items():
      for idx, value2 in value.items():
        for val, cells in value2.items():
          if len(cells) == 1:
            update_flg = True
            cell = cells.pop()
            board.update_cell(cell, val)
    return update_flg

  # ロジック3(仮称)
  # 縦・横・矩形を見て、入る数のリストと入る数が全て同じならば、それ以外のその数リストは候補リストから外す。

  # 力技で解く
  # 全パターンから深さ優先探索で正解を探索する
  # TODO あんまりうまくいっていない。。。探索方法を再検討したい
  def depth_first_search(self, board):
    search_result = []  # 探索結果(値を入れたセルを順に格納するリスト)
    search_history = {}  # keys : セル位置, values : 探索済みセルリスト, その時点での候補リストを格納するタプル
    now_candidate_list = copy.deepcopy(board.candidate_list)

    start_time = time.time()
    # 初期値
    prev_pos = (-1, -1)
    search_history[prev_pos] = ([], copy.deepcopy(board.candidate_list))
    cnt_verbose = 0
    while True: # 完了するまで繰り返す
      cnt_verbose = cnt_verbose + 1
      if cnt_verbose % 100 == 0:
        print("{:10d} ({:2d} / {:2d}) search_result ... {}".format(
          cnt_verbose, len(search_result), board.data._values.flatten().tolist().count(0) + len(search_result), search_result))
      empty_cell_list = board.get_empty_cell_list()
      # print("empty_cell_list : {}".format(empty_cell_list))
      np.random.shuffle(empty_cell_list)
      pos = self.get_next_search_cell_randomly(empty_cell_list, search_history, prev_pos)
      # if board.is_wrong():
      #   print("WRONG...why??? 03")
      if len(pos) == 0: # 次に探索するセルが見つからない場合、探索結果を1つ戻す
        # 次のセルが見つからないため、最新の探索結果のセルを破棄する
        # print("{} ({}, {}) / search_result (minus) ... {}".format(
        #   board.data._values.flatten().tolist().count(0) + len(search_result),
        #   board.data._values.flatten().tolist().count(0), len(search_result), search_result))
        pos_ = search_result.pop(-1)
        # セルの値を0に戻す
        board.update_cell(pos_, 0)
        # 探索履歴を削除する
        search_history.pop(pos_)
        # 1つ前の結果を取り出す
        prev_pos = copy.deepcopy(search_result[-1])
        # 候補リストを元に戻す
        board.candidate_list = copy.deepcopy(search_history[prev_pos][1])
        # print("{} ... {}".format(prev_pos, search_history[prev_pos][0]))
      else:
        # 探索履歴に追加する
        search_history[prev_pos][0].append(pos)
        # 選んだセルに値を入れて検証し、探索結果を更新する
        # if board.is_wrong():
        #   print("WRONG...why??? 01")
        update_flg = self.verify_update_board(board, pos)

        if update_flg: # 更新された
          # if board.is_wrong():
          #   print("WRONG...why??? 02")
          prev_pos = copy.deepcopy(pos)
          search_result.append(pos)
          search_history[pos] = ([], copy.deepcopy(board.candidate_list))
          # print("{} ({}, {}) / search_result (plus) ... {}".format(
          #   board.data._values.flatten().tolist().count(0) + len(search_result),
          #   board.data._values.flatten().tolist().count(0), len(search_result), search_result))

      if board.is_complete(): # 完了した
        end_time = time.time()
        print("Elapsed time[sec] : {:3f}".format(end_time - start_time))
        break

  # ランダムに空いているセルを選択する
  # TODO ランダムに選ぶのが微妙な気がする
  def get_next_search_cell_randomly(self, empty_cell_list, search_history, prev_pos):
    for pos in empty_cell_list:
      if pos not in search_history[prev_pos][0]: # すでに探索済みのセルでない
        return pos
    return ()

  # 探索結果を更新する
  def verify_update_board(self, board, pos):
    # セルに値を入れて検証する
    num_list = copy.deepcopy(board.candidate_list[Const.CELL_KEY_NAME][pos])
    np.random.shuffle(num_list)
    for val in num_list:
      board.update_cell(pos, val)
      prev_candidate_list = copy.deepcopy(board.candidate_list)
      # 候補リストを更新
      board.update_candidate_list(pos, val)

      # 検証
      if board.is_wrong(): # 間違い
        board.update_cell(pos, 0)
        board.candidate_list = copy.deepcopy(prev_candidate_list)
      else: # 間違いでない
        return True

    return False

  # 一番候補が少ない探索するセルを求める。
  def select_search_cell(self, board):
    candidate_list_cell = board.candidate_list[Const.CELL_KEY_NAME]
    min_len_nums = board.size
    # 候補がなるべく少ないセルを選びたい
    cell_list = []
    for cell, nums in candidate_list_cell.items():
      len_nums = len(nums)
      if len_nums == 0:
        continue
      if min_len_nums > len_nums:
        min_len_nums = len_nums
        cell_list.clear()
        cell_list.append(cell)
      elif min_len_nums == len_nums:
        cell_list.append(cell)
    # np.random.shuffle(cell_list)
    # TODO もう少し細かく探索セルを決めたいかも。
    # 候補が少ないセルのうち、他の候補から決定できる個数が多いものを選びたい
    cell = cell_list[np.random.choice(len(cell_list))]
    val = candidate_list_cell[cell][np.random.choice(min_len_nums)]
    return cell, val

max_stars_repo_name: alco-technics/sudoku-solver

max_issues_repo_name: alco-technics/sudoku-solver

max_forks_repo_name: alco-technics/sudoku-solver


text:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on December 2020
@author: Thomas Bonald <bonald@enst.fr>
"""
import unittest

import numpy as np
from numpy.linalg import norm

from sknetwork.data import karate_club
from sknetwork.utils import get_neighbors

class TestNeighbors(unittest.TestCase):

    def test_graph(self):
        adjacency = karate_club()
        neighbors = get_neighbors(adjacency, 5)
        neighbors_true = np.array([0, 6, 10, 16])
        self.assertEqual(norm(neighbors - neighbors_true), 0)

max_stars_repo_name: valeeraZ/scikit-network

max_issues_repo_name: valeeraZ/scikit-network

max_forks_repo_name: valeeraZ/scikit-network

reason: import numpy,from numpy


text:
"""
  Hanging weights puzzle in cpmpy.

  From 'Using LINQ to solve puzzles'
  http://blogs.msdn.com/lukeh/archive/2007/03/19/using-linq-to-solve-puzzles.aspx
  '''
  Here's a puzzle similar to the one in the puzzle hunt.  The diagram
  below is a bunch of weights (A-M) hanging from a system of bars.
  Each weight has an integer value between 1 and 13, and the goal is
  to figure out what each weight must be for the the diagram below to
  balance correctly as shown:

                           |
                           |
               +--+--+--+--+--+--+--+
               |                    |
               |                    |
            +--+--+--+--+--+        |
            |     L        M        |
            |                       |
   +--+--+--+--+--+--+     +--+--+--+--+--+
   H              |  I     |  J        K  |
                  |        |              |
         +--+--+--+--+--+  |     +--+--+--+--+--+
         E              F  |     G              |
                           |                    |
               +--+--+--+--+--+  +--+--+--+--+--+--+
               A              B  C                 D

  The rules for this kind of puzzle are:
  (1) The weights on either side of a given pivot point must be equal,
      when weighted by the distance from the pivot, and
  (2) a bar hanging beneath another contributes it's total weight as
      through it were a single weight.  For instance, the bar on the bottom
      right must have 5*C=D, and the one above it must have 3*G=2*(C+D).
  '''

  Here is a solution

                            |
                            |
                +--+--+--+--+--+--+--+
                |                    |
                |                    |
             +--+--+--+--+--+        |
             |     7        5        |
             |                       |
    +--+--+--+--+--+--+     +--+--+--+--+--+
   11              |  1     |  4        13 |
                   |        |              |
          +--+--+--+--+--+  |     +--+--+--+--+--+
          6              9  |     8              |
                            |                    |
                +--+--+--+--+--+  +--+--+--+--+--+--+
                3              12 2                 10

Model created by Hakan Kjellerstrand, hakank@hakank.com
See also my cpmpy page: http://www.hakank.org/cpmpy/

"""
from cpmpy import *
from cpmpy.solvers import *
import numpy as np
from cpmpy_hakank import *

def show(v):
    return str(v.value())

def hanging_weights():

    N = 13
    x = intvar(1,N,shape=N,name="x")
    a,b,c,d,e,f,g,h,i,j,k,l,m = x
    # total = intvar(1,n*n,name="total")

    model = Model([# total == sum(x),
                   AllDifferent(x),
                   4 * a == b,
                   5 * c == d,
                   3 * e == 2 * f,
                   3 * g == 2 * (c + d),
                   3 * (a + b) + 2 * j == k + 2 * (g + c + d),
                   3 * h == 2 * (e + f) + 3 * i,
                   (h + i + e + f) == l + 4 * m,
                   4 * (l + m + h + i + e + f) == 3 * (j + k + g + a + b + c + d)
                   ])

    ss = CPM_ortools(model)
    # ss.ort_solver.parameters.num_search_workers = 8 # Don't work together with SearchForAllSolutions
    # ss.ort_solver.parameters.search_branching = ort.PORTFOLIO_SEARCH
    # ss.ort_solver.parameters.cp_model_presolve = False
    ss.ort_solver.parameters.linearization_level = 0
    ss.ort_solver.parameters.cp_model_probing_level = 0

    while ss.solve() is not False:
        print(x.value())
        print()
        # print(total.value())
        # This is copied from the MiniZinc model.
        print(
"                          |" + "\n" +
"                          |"+ "\n" +
"              +--+--+--+--+--+--+--+" + "\n" +
"              |                    |" + "\n" +
"              |                    |"+ "\n" +
"           +--+--+--+--+--+        |"+ "\n" +
"           |     "+show(l)+"        " +show(m)+"        |"+ "\n" +
"           |                       |"+ "\n" +
"  +--+--+--+--+--+--+     +--+--+--+--+--+"+ "\n" +
" "+show(h)+"              |  "+show(i)+"     |  "+show(j)+"        "+show(k)+" |"+ "\n" +
"                 |        |              |"+ "\n" +
"        +--+--+--+--+--+  |     +--+--+--+--+--+"+ "\n" +
"        "+show(e)+"              "+show(f)+"  |     "+show(g)+"              |"+ "\n" +
"                          |                    |"+ "\n" +
"              +--+--+--+--+--+  +--+--+--+--+--+--+"+ "\n" +
"              "+show(a)+"              "+show(b)+" "+show(c)+"                 "+show(d)+"\n")

        get_different_solution(ss,x)

hanging_weights()

max_stars_repo_name: hakank/hakank

max_issues_repo_name: hakank/hakank

max_forks_repo_name: hakank/hakank


text:
"""
This is intended to be plug-and-play short squeeze analysis.  Generates top picks and plots.

Need to add in earnings dates into the mix.  Stocks with earnings coming up should be looked into more.
"""
import sys

import talib
import pandas as pd
import numpy as np
from tqdm import tqdm
from sklearn.ensemble import RandomForestRegressor
from sklearn.model_selection import ParameterGrid
import matplotlib.pyplot as plt

# custom
sys.path.append('../stock_prediction/code')
import data_processing as dp

def load_stocks_calculate_short_corr(only_latest=False):
    """
    only_latest: boolean; if True, will keep only the latest 50 days to speed up procesing
    """
    dfs, sh_int, fin_sh = dp.load_stocks(stocks=None)
    if only_latest:
        earliest = sh_int.index[-50]
        sh_int = sh_int[sh_int.index >= earliest]

    latest_stocks = []
    all_sh_stocks = []
    all_sh_stocks_full = []
    latest_date = sh_int['SPY'].index[-1]
    # TODO: make parallel
    for s in tqdm(sh_int.keys()):
        if sh_int[s].shape[0] == 0:
            print(s, 'is empty')
            continue

        if latest_date != sh_int[s].index[-1]:
            print(s, 'is old')
            continue

        if verbose:
            print(s)

        df = sh_int[s].copy()
        # create future price changes
        for i in [5, 10, 20, 40]:
            df['{}d_future_price_change'.format(i)] = df['Adj_Close'].pct_change(i).shift(-i)

        df['ticker'] = s

        # create short-close correlations -- TODO: deal with -1s
        # if short % is all -1 or 0, won't work.  if less than 20 samples, rolling corr with 20 period window won't work
        # also broke on DF with 22 samples
        if df['Short_%_of_Float'].mean() in [-1, 0] or df.shape[0] < 30:
            df['Short_%_of_Float_10d_EMA'] = -np.inf
            df['Adj_Close_10d_EMA'] = talib.EMA(df['Adj_Close'].values, timeperiod=10)
            df['short_close_corr_10d_EMA'] = -np.inf
            df['short_close_corr_rocr_20d'] = -np.inf
            df['short_%_rocr_20d'] = -np.inf
        else:
            df['Short_%_of_Float_10d_EMA'] = talib.EMA(df['Short_%_of_Float'].values, timeperiod=10)
            df['Adj_Close_10d_EMA'] = talib.EMA(df['Adj_Close'].values, timeperiod=10)

            # essentially, we want to take an arbitrary number of points, calculate correlation, and find where the correlations are largest
            # take 10 points at a time and get correlations first, then take parts that have largest correlations, and keep expanding by 5 points at a time until correlation decreases
            corr = df[['Short_%_of_Float_10d_EMA', 'Adj_Close_10d_EMA']].rolling(window=20).corr()
            df['short_close_corr_10d_EMA'] = corr.unstack(1)['Short_%_of_Float_10d_EMA']['Adj_Close_10d_EMA']
            df['short_close_corr_10d_EMA'].replace(np.inf, 1, inplace=True)
            df['short_close_corr_10d_EMA'].replace(-np.inf, -1, inplace=True)
            df['short_close_corr_10d_EMA'].clip(lower=-1, upper=1, inplace=True)

            # WARNING: things with small (< 1%) Short % of float will result in huge rocr...maybe do something about this
            df['short_close_corr_rocr_20d'] = talib.ROCR100(df['short_close_corr_10d_EMA'].values, timeperiod=20)
            df['short_%_rocr_20d'] = talib.ROCR100(df['Short_%_of_Float_10d_EMA'].values, timeperiod=20)

            # auto-detect long stretches of negative and positive correlation
            thresh = 0.7
            rolling = df['short_close_corr_10d_EMA'].rolling(window=20).min()
            df['Short_%_positive_corr_detection'] = rolling > thresh
            df['Short_%_positive_corr_detection'] = df['Short_%_positive_corr_detection'].astype('int16')
            # sh_int[ticker]['Short_%_positive_corr_detection'].plot()
            # plt.show()

            df['Short_%_negative_corr_detection'] = rolling < -thresh
            df['Short_%_negative_corr_detection'] = df['Short_%_negative_corr_detection'].astype('int16')

        latest_stocks.append(df.iloc[-1])
        all_sh_stocks_full.append(df)
        all_sh_stocks.append(df.dropna())

    latest_stocks_df = pd.concat(latest_stocks, axis=1).T
    latest_stocks_df.set_index('ticker', inplace=True)

    all_sh_stocks_df = pd.concat(all_sh_stocks)
    all_sh_stocks_df['market_cap'] = all_sh_stocks_df['Shares_Outstanding'] * all_sh_stocks_df['Adj_Close']
    all_sh_stocks_full_df = pd.concat(all_sh_stocks_full)

    return all_sh_stocks_df, all_sh_stocks_full_df, latest_stocks_df, sh_int

def make_larger_shorts(all_sh_stocks_df):
    """
    gets dataframe of stocks with larger amounts of shares short overall (not only currently)
    """
    # plot short close rocr 20d to see if trend from TDOC holds overall -- under 80% ROCR100 means high gains
    # no overall trend, but for individual stocks there tends to be a trend
    shorts = all_sh_stocks_df[all_sh_stocks_df['short_%_rocr_20d'] != -np.inf]
    shorts_nona = shorts.dropna()

    ticker_groups = shorts_nona[['ticker', 'Short_%_of_Float', 'Days_to_Cover', 'rocp_cl', 'short_%_rocr_20d', 'short_close_corr_rocr_20d']].groupby('ticker').mean()
    larger_shorts = ticker_groups[ticker_groups['Short_%_of_Float'] > 10]
    short_stocks = larger_shorts.index

    return ticker_groups, larger_shorts, short_stocks

def short_squeeze_analysis(all_sh_stocks_df, ticker='TDOC'):
    # also look at EMA
    df = all_sh_stocks_df[all_sh_stocks_df['ticker'] == ticker]

    # I think all this is already calcuclated in load_stocks_calculate_short_corr
    # df['Short_%_of_Float_10d_EMA'] = talib.EMA(sh_int[ticker]['Short_%_of_Float'].values, timeperiod=10)
    # df['Adj_Close_10d_EMA'] = talib.EMA(df['Adj_Close'].values, timeperiod=10)
    df[['Short_%_of_Float', 'Short_%_of_Float_10d_EMA']].plot()

    # look for sections where short % mva and price are (ideally, inversely) correlated
    # make gradient of color where oldest -> newest goes from dark to light
    cm = plt.get_cmap('copper')
    num_colors = df.shape[0]
    colors = []
    for i in range(num_colors):
        colors.append(cm(i / num_colors))  # color will now be an RGBA tuple

    df.plot.scatter(x='Short_%_of_Float_10d_EMA', y='Adj_Close', color=colors)
    plt.show()

    # essentially, we want to take an arbitrary number of points, calculate correlation, and find where the correlations are largest
    # take 10 points at a time and get correlations first, then take parts that have largest correlations, and keep expanding by 5 points at a time until correlation decreases
    # already being done in preprocessing now
    # corr = df[['Short_%_of_Float_10d_EMA', 'Adj_Close_10d_EMA']].rolling(window=20).corr()
    # df['short_close_corr_10d_EMA'] = corr.unstack(1)['Short_%_of_Float_10d_EMA']['Adj_Close_10d_EMA']
    # df['short_close_corr_10d_EMA'].replace(np.inf, 1, inplace=True)
    # df['short_close_corr_10d_EMA'].replace(-np.inf, -1, inplace=True)
    # df['short_close_corr_10d_EMA'].clip(lower=-1, upper=1, inplace=True)

    # rows, columns
    # f, ax = plt.subplots(nrows=2, ncols=1, sharex=True)
    # sh_int[ticker]['short_close_corr_10d_EMA'].plot(ax=ax[0])
    # sh_int[ticker][['Adj_Close_10d_EMA', 'Short_%_of_Float_10d_EMA']].plot(ax=ax[1])
    # plt.show()

    # auto-detect long stretches of negative and positive correlation
    # now calculated in load_stocks_calculate_short_corr function
    # thresh = 0.8
    # rolling = df['short_close_corr_10d_EMA'].rolling(window=20).min()
    # df['Short_%_positive_corr_detection'] = rolling > thresh
    # df['Short_%_positive_corr_detection'] = df['Short_%_positive_corr_detection'].astype('int16')
    # # sh_int[ticker]['Short_%_positive_corr_detection'].plot()
    # # plt.show()
    #
    # df['Short_%_negative_corr_detection'] = rolling < -thresh
    # df['Short_%_negative_corr_detection'] = df['Short_%_negative_corr_detection'].astype('int16')
    df[['short_close_corr_10d_EMA', 'Short_%_negative_corr_detection', 'Adj_Close', 'Short_%_of_Float_10d_EMA', 'Days_to_Cover', 'short_close_corr_rocr_20d']].plot(subplots=True)
    plt.show()

all_sh_stocks_df, all_sh_stocks_full_df, latest_stocks_df, sh_int = load_stocks_calculate_short_corr()
ticker_groups, larger_shorts, short_stocks = make_larger_shorts(all_sh_stocks_df)

short_squeeze_analysis(all_sh_stocks_full_df, ticker='TLRY')


text:
"""
Script for making init generation codes from given images
"""

import os
import numpy as np
import Generators
from caffe.io import load_image
from cv2 import imwrite
from shutil import copyfile

# parameters
generator_name = 'deepsim-conv4'
generator = Generators.get_generator(generator_name)
srcdirs = ('imgs',)    # to be changed; list of directories containing source images
dstrootdir = os.path.join('codes', generator_name)

# main
for srcdir in srcdirs:
    leafdir = os.path.basename(srcdir)
    dstdir = os.path.join(dstrootdir, leafdir)
    if not os.path.isdir(dstdir):
        os.makedirs(dstdir)
    imgfns = [fn for fn in os.listdir(srcdir) if len(fn) > 3 and fn[-4:] == '.bmp']
    for imgfn in imgfns:
        im = load_image(os.path.join(srcdir, imgfn))
        if generator.name == 'raw_pixel':
            code = generator.encode(im)
        else:
            code = generator.encode(im, steps=200)
        np.save(os.path.join(dstdir, imgfn[:-4]+'.npy'), code, allow_pickle=False)
        imwrite(os.path.join(dstdir, imgfn[:-4]+'_syn.png'), generator.visualize(code)[:, :, ::-1])
        copyfile(os.path.join(srcdir, imgfn), os.path.join(dstdir, imgfn))

max_stars_repo_name: willwx/XDream

max_issues_repo_name: willwx/XDream

max_forks_repo_name: willwx/XDream


text:
import os
import numpy as np
from keras.utils import np_utils
import models
import extract_feats.opensmile as of
import extract_feats.librosa as lf
import utils.opts as opts

'''
train(): 训练模型

输入:
	config(Class)

输出:
    model: 训练好的模型
'''
def train(config):

    # 加载被 preprocess.py 预处理好的特征
    if(config.feature_method == 'o'):
        x_train, x_test, y_train, y_test = of.load_feature(config, config.train_feature_path_opensmile, train = True)

    elif(config.feature_method == 'l'):
        x_train, x_test, y_train, y_test = lf.load_feature(config, config.train_feature_path_librosa, train = True)

    # x_train, x_test (n_samples, n_feats)
    # y_train, y_test (n_samples)

    # 搭建模型
    model = models.setup(config = config, n_feats = x_train.shape[1])

    # 训练模型
    print('----- start training', config.model, '-----')
    if config.model in ['lstm', 'cnn1d', 'cnn2d']:
        y_train, y_val = np_utils.to_categorical(y_train), np_utils.to_categorical(y_test) # 独热编码
        model.train(
            x_train, y_train,
            x_test, y_val,
            batch_size = config.batch_size,
            n_epochs = config.epochs
        )
    else:
        model.train(x_train, y_train)
    print('----- end training ', config.model, ' -----')

    # 验证模型
    model.evaluate(x_test, y_test)
    # 保存训练好的模型
    model.save_model(config)

if __name__ == '__main__':

    config = opts.parse_opt()
    train(config)


text:
"""
@File : recognition.py
@Author: Dong Wang
@Date : 2020/5/5
"""
from . import load_model, gestures, Sub2Model
from PIL import Image
import numpy as np
import six
import os

class recognition(object):
    def __init__(self, model):
        """
        Constructor of the abstract class. It is designed for common interface.
        You should implement a subclass of recognition firstly.

        # Usage:

        model = load_model(modelPath) # modelPath is a trained neural model.
        gR = gestureRecognition(model)
        for landmark in landmarks:
            gesture = gR.recognize(landmark)
            print("The gesture is ", gesture)

        """
        self.model = model

    def recognize(self, landmark):
        """ recognize gesture from an image
        :param landmark: landmark array
        :return:gesture
        """
        # landmark = landmark
        # # Do something to the landmark
        # gesture = self.model.predict(landmark)
        # return gesture
        raise NotImplementedError

class gestureRecognition(recognition):
    def __init__(self, model=Sub2Model):
        model = load_model(model)
        super(gestureRecognition, self).__init__(model=model)

    def recognize(self, landmark):
        response = {'code': 200}
        if len(landmark) != 80:
            response['code'] = 401
            return response
        elif isinstance(landmark, np.ndarray):
            # Img is Image type
            response['code'] = 200
            response['gesture'], response['probability'] = self.predict(landmark)
        else:
            response['code'] = 402
        return response

    def predict(self, pointsArray):
        if len(pointsArray.shape) == 1:
            pointsArray = np.array([pointsArray])
        pred = self.model.predict(pointsArray)[0]
        index = int(np.argmax(pred))
        gesture = gestures[index]
        return gesture, pred[index]


text:
from io import StringIO
import numpy as np
import inspect
import yaml

def code_print(code: str, file=None, print_function=None):
    """
    Print function used print out code with line numbers.

    :param code: String containing the code to print.
    :param file:
    :param print_function: Print function to use. By default this
                           is the python print(...) function.
    :return:
    """
    str_io = StringIO()

    lines = code.splitlines()
    n_chars = int(np.floor(np.log10(len(lines))+1))

    for i, line in enumerate(lines):
        str_io.write(f"{i+1:{n_chars}}| {line}\n")

    (print_function or print)(str_io.getvalue(), file=file)

def debug_print(*args, sep=' ', file=None):
    """
    Print function that appends the printed text with the location of the print statement.

    :param args:
    :param sep:
    :param file:
    :return:
    """
    file_name = inspect.stack()[1][1]
    line_no = inspect.stack()[1][2]
    function_name = inspect.stack()[1][3]
    msg = sep.join(map(str, args))
    print(f'{msg}: File "{file_name}", line {line_no}, in {function_name}', file=file)

def verbose_print(verbose, print_function=None):
    """
    Factory function that return a working print function if verbose is True.
    Otherwise it returns a dummy function that won't print anything.

    :param verbose: True to return print_function
    :param print_function: Optional print function to use. Uses 'print' by default.
    :return:
    """

    if verbose:
        return print_function or print
    else:
        def vprint(*args, **kwars):
            pass
        return vprint

def yaml_print(obj, default_flow_style=False, print_func=debug_print):
    """

    :param obj: Object to print
    :param default_flow_style: True to use PyYAML's default flow style. For more
                               information check PyYAML's documentation.
    :param print_func: Print function to use. By default its debug print.
    :return:
    """
    return print_func(yaml.dump(obj, default_flow_style=default_flow_style))

# Monkey patch debug print and yaml_print into the buildins so they are available everywhere.
__builtins__["debug_print"] = debug_print
__builtins__["yaml_print"] = yaml_print

max_stars_repo_name: aarondewindt/cw

max_issues_repo_name: aarondewindt/cw

max_forks_repo_name: aarondewindt/cw


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

from __future__ import division
from math import pi, cos, sin, atan, acos
import numpy as np

class Circle(object):
    def __init__(self, p, r):
        self.c = p
        self.r = float(r)

    def __repr__(self):
        return self.__str__()

    def __str__(self):
        return "".join(['Circle[' + self.c.__str__() + ':' + str(self.r) + ']'])

    def __eq__(self, other):
        return self.c == other.c and self.r == other.r

class Point(object):
    def __init__(self, x, y=0, z=0):
        if(type(x) in [list, tuple, np.ndarray]):
            self.x = x[0] if len(x) >= 1 else 0.0
            self.y = x[1] if len(x) >= 2 else 0.0
            self.z = x[2] if len(x) >= 3 else 0.0
        else:
            self.x = x
            self.y = y
            self.z = z

    def __str__(self):
        return "".join(['Point(' + str(self.x) + ',' + str(self.y) + ',' + str(self.z) + ')'])

    def __eq__(self, other):
        return self.__dict__ == other.__dict__

    def __sub__(self, other):
        if isinstance(other, Point):
            tx = self.x - other.x
            ty = self.y - other.y
            tz = self.z - other.z
        else:
            tx = self.x - other.dx
            ty = self.y - other.dy
            tz = self.z - other.dz
        return Point(tx, ty, tz)

    def __add__(self, other):
        if isinstance(other, Point):
            tx = self.x + other.x
            ty = self.y + other.y
            tz = self.z + other.z
        else:
            tx = self.x + other.dx
            ty = self.y + other.dy
            tz = self.z + other.dz
        return Point(tx, ty, tz)

    def __mul__(self, other):
        return Point(other * self.x, other * self.y, other * self.z)

    def __rmul__(self, other):
        return Point(other * self.x, other * self.y, other * self.z)

    def __div__(self, other):
        return Point(self.x / other, self.y / other, self.z / other)

    def __neg__(self):
        x = -self.x
        y = -self.y
        z = -self.z
        return Point(x, y, z)

    def __getitem__(self, id):
        return [self.x, self.y, self.z][id]

    def area(self):
        return 0.0

    def dist(self, other):
        return ((self[0] - other[0]) ** 2 + (self[1] - other[1]) ** 2 + (self[2] - other[2]) ** 2) ** 0.5

    def std(self):
        return [self.x, self.y, self.z]

    def c2s(self):
        R = self.dist(Point(0, 0, 0))
        lg = atan(self.y / self.x)
        lat = acos(self.z / R)
        return (lg, lat, R)

    def transform(self, p, rot):
        px = cos(rot) * self.x + sin(rot) * self.y
        py = -sin(rot) * self.x + cos(rot) * self.y
        p_t = Point(px, py)
        return p_t - p

    def rot(self, a):
        px = cos(a) * self.x - sin(a) * self.y
        py = sin(a) * self.x + cos(a) * self.y
        return Point(px, py)

max_stars_repo_name: AlexisTM/MultilaterationTDOA

max_forks_repo_name: AlexisTM/MultilaterationTDOA