text:
# The Leginon software is Copyright 2004
# The Scripps Research Institute, La Jolla, CA
# For terms of the license agreement
# see http://ami.scripps.edu/software/leginon-license
#
# $Source: /ami/sw/cvsroot/pyleginon/navigator.py,v $
# $Revision: 1.130 $
# $Name: not supported by cvs2svn $
# $Date: 2008-02-14 22:08:15 $
# $Author: pulokas $
# $State: Exp $
# $Locker:  $

# leginon
import node
import event
from leginon import leginondata
import calibrationclient
import gui.wx.Navigator
import instrument
import presets
import cameraclient

# myami
from pyami import correlator, peakfinder, imagefun, ordereddict
import pyami.fileutil
import pyami.mrc

# python standard
import threading
import time
import types
import os.path
import itertools
import math
import logging

# 3rd party
import numpy
import scipy.ndimage

# used to create debug filenames
idcounter = itertools.cycle(range(100))

class NavigatorClient(object):
	eventoutputs = [event.MoveToTargetEvent]
	eventinputs = [event.MoveToTargetDoneEvent]

	def __init__(self, node):
		self.node = node
		self.node.addEventInput(event.MoveToTargetDoneEvent, self.handleMoveDone)
		self.movedone = threading.Event()

	def handleMoveDone(self, evt):
		self.movedonestatus = evt['status']
		self.movedone.set()

	def moveToTarget(self, target, movetype, precision=0.0, accept_precision=1e-3, final_imageshift=False, use_current_z=False):
		'''
		Go through a thread of moveToTarget in Navigator
		'''
		self.node.startTimer('moveToTarget')
		ev = event.MoveToTargetEvent(target=target, movetype=movetype)
		ev['move precision'] = precision
		ev['accept precision'] = accept_precision
		ev['final image shift'] = final_imageshift
		ev['use target z'] = not use_current_z
		self.movedone.clear()
		self.node.outputEvent(ev, wait=False)
		## wait for event
		self.movedone.wait()
		self.node.stopTimer('moveToTarget')
		return self.movedonestatus

class Navigator(node.Node):
	panelclass = gui.wx.Navigator.Panel
	settingsclass = leginondata.NavigatorSettingsData
	defaultsettings = {
		'instruments': {'tem':None, 'ccdcamera':None},
		'pause time': 2.5,
		'move type': 'image shift',
		'check calibration': True,
		'precision': 0.0,
		'accept precision': 1e-3,
		'complete state': True,
		'override preset': False,
		'max error': 256,
		'cycle after': False,
		'cycle each': False,
		'final image shift': False,
		'camera settings': cameraclient.default_settings,
		'preexpose': True,
	}
	eventinputs = node.Node.eventinputs + presets.PresetsClient.eventinputs + [event.MoveToTargetEvent]
	eventoutputs = node.Node.eventoutputs + presets.PresetsClient.eventoutputs + [event.CameraImagePublishEvent, event.MoveToTargetDoneEvent,event.UpdatePresetEvent]

	def __init__(self, id, session, managerlocation, **kwargs):
		node.Node.__init__(self, id, session, managerlocation, **kwargs)
		self.instrument = instrument.Proxy(self.objectservice, self.session,
																				self.panel)
		self.calclients = ordereddict.OrderedDict()
		self.calclients['image shift'] = calibrationclient.ImageShiftCalibrationClient(self)
		self.calclients['stage position'] = calibrationclient.StageCalibrationClient(self)
		self.calclients['modeled stage position'] = calibrationclient.ModeledStageCalibrationClient(self)
		self.calclients['beam shift'] = calibrationclient.BeamShiftCalibrationClient(self)
		self.calclients['image beam shift'] = calibrationclient.ImageBeamShiftCalibrationClient(self)

		self.pcal = calibrationclient.PixelSizeCalibrationClient(self)
		self.presetsclient = presets.PresetsClient(self)
		self.stagelocations = []
		self.getLocationsFromDB()

		self.correlator = correlator.Correlator()
		self.peakfinder = peakfinder.PeakFinder()
		self.oldimagedata = None
		self.newimagedata = None
		self.oldstate = None
		self.newstate = None

		self.addEventInput(event.MoveToTargetEvent, self.handleMoveEvent)

		self.start()

	def handleMoveEvent(self, ev):
		nodename = ev['node']
		movetype = ev['movetype']
		targetdata = ev['target']
		precision = ev['move precision']
		accept_precision = ev['accept precision']
		final_imageshift = ev['final image shift']
		use_target_z = ev['use target z']
		self.logger.info('handling %s request from %s' % (movetype, nodename,))
		imagedata = targetdata['image']
		self.newImage(imagedata)
		preset = imagedata['preset']
		rows = targetdata['delta row']
		cols = targetdata['delta column']

		if precision:
			check=True
		else:
			check=False

		## initial change to parent preset,
		## otherwise presets manager may never know what happened
		self.logger.info('change to parent preset: %s' % (preset['name'],))
		self.presetsclient.toScope(preset['name'])

		self.startTimer('move')
		# Force cycle_after to True because PresetsManager does not know that preset
		# has been changed by Navigator and will not cycle on the first target.
		# Later targets in the list do cycle according to PresetsManager settings
		# even if Navigator settings['cycle after'] is off
		if use_target_z:
			# if the move comes from target adjustment, z focus is likely done.
			# Therefore target z should not be set according to parent
			msg = 'Set to target z %.6f' % imagedata['scope']['stage position']['z']
			self.testprint('Navigator: ' + msg)
			self.logger.debug(msg)
			self.instrument.tem.setStagePosition({'z':imagedata['scope']['stage position']['z']})
		stagenow = self.instrument.tem.StagePosition
		self.logger.debug('Navigator: z in navigator move %.6f' % stagenow['z'])
		status = self.move(rows, cols, movetype, precision, accept_precision, check, preset=preset, final_imageshift=final_imageshift, cycle_after=True)
		self.stopTimer('move')

		evt = event.MoveToTargetDoneEvent(status=status, target=targetdata)
		evt['destination'] = nodename
		self.outputEvent(evt)

	def newImage(self, imagedata):
		self.oldimagedata = self.newimagedata
		self.newimagedata = imagedata
		self.setImage(imagedata['image'], 'Image')
		self.correlator.insertImage(imagedata['image'])

	## called by GUI when image is clicked
	def navigate(self, xy):
		movetype = self.settings['move type']
		precision = self.settings['precision']
		accept_precision = self.settings['accept precision']
		final_imageshift = self.settings['final image shift']
		clickrow = xy[1]
		clickcol = xy[0]
		clickshape = self.newimagedata['image'].shape

		# calculate delta from image center
		centerr = clickshape[0] / 2.0 - 0.5
		centerc = clickshape[1] / 2.0 - 0.5
		deltarow = clickrow - centerr
		deltacol = clickcol - centerc

		check = self.settings['check calibration']
		status = self.move(deltarow, deltacol, movetype, precision, accept_precision, check, final_imageshift=final_imageshift)

		## acquire image if check not done
		if not check:
			self.reacquireImage()
		elif status != 'error' and movetype != 'image shift' and final_imageshift:
			self.reacquireImage()

		self.panel.navigateDone()

	def _moveback(self):
		self.logger.info('moving back to previous state')
		emdat = leginondata.NavigatorScopeEMData()
		emdat.friendly_update(self.oldstate)
		try:
			self.instrument.setData(emdat)
		except:
			self.logger.exception('unable to set instrument')
			return

	def _move(self, row, col, movetype):
		errstr = 'Move failed: %s'
		# get relavent info from click event

		self.logger.info('Moving...')

		pixelshift = {'row':-row, 'col':-col}
		scope = self.newimagedata['scope']
		camera = self.newimagedata['camera']

		# figure out shift
		calclient = self.calclients[movetype]
		try:
			newstate = calclient.transform(pixelshift, scope, camera)
		except calibrationclient.NoMatrixCalibrationError, e:
			errsubstr = 'unable to find calibration for %s' % e
			self.logger.error(errstr % errsubstr)
			self.beep()
			self.panel.navigateDone()
			return True
		except Exception, e:
			self.logger.exception(errstr % e)
			self.beep()
			self.panel.navigateDone()
			return True

		# Display the movement applied
		if movetype.endswith('stage position'):
			moveparams = ['stage position',]
		elif movetype == 'image beam shift':
			moveparams = ['image shift','beam shift']
		else:
			moveparams = [movetype,]
		for moveparam in moveparams:
			scopeshift = {'moveparam': moveparam, 'x':None, 'y':None}
			for axis in ('x','y'):
				scopeshift[axis] = newstate[moveparam][axis] - scope[moveparam][axis]
			self.logger.info('change in %(moveparam)s: %(x).4e, %(y).4e' % scopeshift)

		# Avoid changing stage z according to imagedata here since it is handled
		# before the move
		stagenow = self.instrument.tem.StagePosition
		newstate['stage position']['z'] = stagenow['z']

		self.oldstate = self.newstate
		self.newstate = newstate
		emdat = leginondata.NavigatorScopeEMData()
		emdat.friendly_update(newstate)
		try:
			self.instrument.setData(emdat)
		except:
			self.logger.exception(errstr % 'unable to set instrument')
			return True

		return False

	def move_away_move_back(self, label, moves, distance, angle=None):
		'''
		Test reproducibility of the stage from position of a given distance.
		'''
		self.logger.info('start reproducibility test')
		# initial move
		origstage = self.instrument.tem.StagePosition
		origx = origstage['x']
		origy = origstage['y']
		if angle is None:
			rangle = 2 * numpy.pi * numpy.random.rand()
		else:
			rangle = angle
		deltax = distance * numpy.cos(rangle)
		deltay = distance * numpy.sin(rangle)
		tmpx = origx + deltax
		tmpy = origy + deltay
		self.instrument.tem.StagePosition = {'x': tmpx, 'y': tmpy}
		time.sleep(2)
		self.instrument.tem.StagePosition = {'x': origx, 'y': origy}
		time.sleep(2)
		self.acquireImage()

		self.origmove = 0,0
		self.origimagedata = self.newimagedata
		origscope = self.origimagedata['scope']
		origcamera = self.origimagedata['camera']
		origstage = origscope['stage position']
		origx = origstage['x']
		origy = origstage['y']
		for i in range(moves):
			if angle is None:
				rangle = 2 * numpy.pi * numpy.random.rand()
			else:
				rangle = angle
			deltax = distance * numpy.cos(rangle)
			deltay = distance * numpy.sin(rangle)
			tmpx = origx + deltax
			tmpy = origy + deltay
			self.instrument.tem.StagePosition = {'x': tmpx, 'y': tmpy}
			time.sleep(2)
			self.instrument.tem.StagePosition = {'x': origx, 'y': origy}
			time.sleep(2)
			self.reacquireImage()
			r,c,dist = self.checkMoveError()
			self.logger.info('move error: pixels: %s, %s, %.3em,' % (r,c,dist,))
			repdata = leginondata.StageReproducibilityData()
			repdata['session'] = self.session
			repdata['label'] = label
			repdata['move x'] = deltax
			repdata['move y'] = deltay
			repdata['error pixels r'] = r
			repdata['error pixels c'] = c
			repdata['error meters'] = dist
			repdata.insert(force=True)
		self.logger.info('reproducibility test done')

	def cycleToPreset(self, preset=None, keep_shift=False):
		if preset is None:
			return
		self.logger.info('preset cycle')
		self.presetsclient.toScope(preset['name'], keep_shift=keep_shift)

	def move(self, row, col, movetype, precision=0.0, accept_precision=1e-3, check=False, preset=None, cycle_after=True, final_imageshift=False):
		self.setStatus('processing')
		self.origimagedata = self.newimagedata
		self.origmove = row,col
		err = self._move(row, col, movetype)
		if err:
			self.setStatus('idle')
			return 'error'

		shape = self.origimagedata['image'].shape
		target = shape[0]/2.0-0.5+self.origmove[0], shape[1]/2.0-0.5+self.origmove[1]
		status = 'ok'
		if check:
			# pre-expose after first move and before reacquire
			if preset and self.settings['preexpose'] and preset['pre exposure']:
				self.exposeSpecimenWithShutterOverride(preset['pre exposure'])

			if self.outofbounds(target, shape):
				self.logger.info('target out of bounds, so cannot check error')
				self.setStatus('idle')
				self.reacquireImage()
				return 'error'
			if self.settings['cycle each']:
				self.cycleToPreset(preset)
			self.reacquireImage()
			r,c,dist = self.checkMoveError()
			self.logger.info('move error: pixels: %s, %s, %.3em,' % (r,c,dist,))

			if precision:
				self.logger.info('checking that move error is less than %.3e' % (precision,))
				while dist > precision:
					time.sleep(0.2)
					self._move(r, c, movetype)
					time.sleep(2.0)
					if self.settings['cycle each']:
						self.cycleToPreset(preset)
					self.reacquireImage()
					lastdist = dist
					r,c,dist = self.checkMoveError()
					self.logger.info('move error: pixels: %s, %s, %.3em,' % (r,c,dist,))
					if dist > lastdist:
						self.logger.info('error got worse')
						if dist > accept_precision:
							status = 'error'
						else:
							self._moveback()
						break

				# final image shift
				if status != 'error' and 'shift' not in movetype and final_imageshift:
					self.logger.info('making final correction with image shift')
					self._move(r, c, 'image shift')

				self.logger.info('move correction done, status: %s' % (status,))
				if self.settings['cycle after'] or cycle_after:
					self.cycleToPreset(preset, keep_shift=True)

		self.setStatus('idle')
		return status

	def outofbounds(self, coord, shape):
		if coord[0] < 0 or coord[0] > shape[0]-1:
			return True
		if coord[1] < 0 or coord[1] > shape[1]-1:
			return True
		return False

	def checkMoveError(self):
		#maxerror = self.settings['max error']
		#limit = (int(maxerror*2), int(maxerror*2))

		oldshape = self.oldimagedata['image'].shape
		limit = oldshape
		location = oldshape[0]/2.0-0.5+self.origmove[0], oldshape[1]/2.0-0.5+self.origmove[1]

		im1 = imagefun.crop_at(self.origimagedata['image'], location, limit, mode='constant', cval=0.0)
		im2 = imagefun.crop_at(self.newimagedata['image'], 'center', limit)

		pc = correlator.phase_correlate(im2, im1, zero=False)
		pc = scipy.ndimage.gaussian_filter(pc,1)
		subpixelpeak = self.peakfinder.subpixelPeak(newimage=pc, guess=(0.5,0.5), limit=limit)
		res = self.peakfinder.getResults()
		unsignedpixelpeak = res['unsigned pixel peak']
		peaktargets = [(unsignedpixelpeak[1], unsignedpixelpeak[0])]
		r_error = subpixelpeak[0]
		c_error = subpixelpeak[1]

		self.setImage(pc, 'Correlation')
		peaktargets = [(unsignedpixelpeak[1], unsignedpixelpeak[0])]
		self.setTargets(peaktargets, 'Peak')

		## calculate error distance
		scope = self.newimagedata['scope']
		camera = self.newimagedata['camera']
		mag = scope['magnification']
		tem = scope['tem']
		ccdcamera = camera['ccdcamera']
		pixelsize = self.pcal.retrievePixelSize(tem, ccdcamera, mag)
		cbin = camera['binning']['x']
		rbin = camera['binning']['y']
		rpix = r_error * rbin
		cpix = c_error * cbin
		pixdist = math.hypot(rpix,cpix)
		distance = pixdist * pixelsize

		return r_error, c_error, distance

	def uiAcquireImage(self):
		self.acquireImage()
		self.panel.acquisitionDone()

	def reacquiredFilename(self):
		orig_id = self.origimagedata.dbid
		if orig_id is None:
			orig_id = 0
		sessionname = self.session['name']
		timestamp = time.strftime('%d%H%M%S', time.localtime())
		nextid = idcounter.next()
		f = '%08d_%s_%02d.mrc' % (orig_id, timestamp, nextid)
		return f

	def reacquireImage(self):
		if self.origimagedata is None:
			raise RuntimeError('no image to reacquire')
		# configure camera just like current image
		cameradata = self.origimagedata['camera']
		ccdcamera = cameradata['ccdcamera']
		## except correction channel should be opposite
		if self.origimagedata['correction channel']:
			corchannel = 0
		else:
			corchannel = 1
		#camerasettings = leginondata.CameraSettingsData(initializer=camera)
		self.instrument.setCCDCamera(ccdcamera['name'])
		self.instrument.setData(cameradata)
		self._acquireImage(corchannel)

		# save the reacquired image
		if self.logger.isEnabledFor(logging.DEBUG):
			debugpath = os.path.join(self.session['image path'], 'debug')
			pyami.fileutil.mkdirs(debugpath)
			filename = self.reacquiredFilename()
			filename = os.path.join(debugpath, filename)
			pyami.mrc.write(self.newimagedata['image'], filename)
			self.logger.debug('saving reacquired mrc: %s' %(filename,))

	def acquireImage(self):
		if self.settings['override preset']:
			## use override
			instruments = self.settings['instruments']
			try:
				self.instrument.setTEM(instruments['tem'])
				self.instrument.setCCDCamera(instruments['ccdcamera'])
			except ValueError, e:
				self.logger.error('Cannot set instruments: %s' % (e,))
				return
			except instrument.NotAvailableError, e:
				self.logger.error('%s' % (e,))
			try:
				self.instrument.ccdcamera.Settings = self.settings['camera settings']
			except Exception, e:
				errstr = 'Acquire image failed: %s'
				self.logger.error(errstr % e)
				return
		else:
			# default to current camera config set by preset
			if self.presetsclient.getCurrentPreset() is None:
				self.logger.error('Preset is unknown and preset override is off')
				return
		self._acquireImage()

	def _restoreSaveFrames(self):
			# Bug #3614 fi
			try:
				self.instrument.ccdcamera.SaveRawFrames = self.was_save_frames
			except:
				pass

	def _setWasSaveFrames(self):
		try:
			self.was_save_frames = self.instrument.ccdcamera.SaveRawFrames
		except:
			self.was_save_frames = False

	def _acquireImage(self, channel=0):
		#Must set to the camera to be used by now.
		self._setWasSaveFrames()
		try:
			self.instrument.ccdcamera.SaveRawFrames = False
		except:
			pass
		time.sleep(self.settings['pause time'])
		try:
			imagedata = self.acquireCorrectedCameraImageData(channel=channel)
		except:
			self.logger.error('unable to get corrected image')
			self._restoreSaveFrames()
			return

		if imagedata is None:
			self.logger.error('Acquire image failed')
			self._restoreSaveFrames()
			return

		self._restoreSaveFrames()
		self.newImage(imagedata)
		return imagedata

	def fromScope(self, name, comment='', xyonly=True):
		'''
		create a new location with name
		if a location by this name already exists in my
		list of managed locations, it will be replaced by the new one
		also returns the new location object
		'''
		errstr = 'Location from instrument failed: %s'
		try:
			allstagedata = self.instrument.tem.StagePosition
		except:
			self.logger.error(errstr % 'unable to get stage position')
			return
		stagedata = {}
		stagedata['x'] = allstagedata['x']
		stagedata['y'] = allstagedata['y']
		if not xyonly:
			stagedata['z'] = allstagedata['z']
			stagedata['a'] = allstagedata['a']

		loc = self.getLocation(name)
		if loc is None:
			newloc = leginondata.StageLocationData()
			newloc['name'] = name
			newloc['xy only'] = xyonly
		else:
			newloc = leginondata.StageLocationData(initializer=loc.toDict())

		newloc['session'] = self.session
		newloc.update(stagedata)

		## save comment, remove the old location
		for loc in self.stagelocations:
			if loc['name'] == name:
				comment = loc['comment']
				self.stagelocations.remove(loc)
				break
		newloc['comment'] = comment
		newloc['removed'] = False

		self.stagelocations.append(newloc)

		self.locationToDB(newloc)
		locnames = self.locationNames()
		self.panel.locationsEvent(locnames)
		#self.logger.info('Location names %s' % (locnames,))
		return newloc

	def removeLocation(self, loc):
		'''
		remove a location by index or reference
		loc is either a location, or index of the location
		'''
		locremove = None
		if type(loc) is int:
			locremove = self.stagelocations[loc]
			del self.stagelocations[loc]
		elif type(loc) in types.StringTypes:
			loccopy = list(self.stagelocations)
			for location in loccopy:
				if location['name'] == loc:
					locremove = location
					self.stagelocations.remove(location)
		else:
			locremove = loc
			self.stagelocations.remove(loc)

		if locremove is not None:
			removeloc = leginondata.StageLocationData(initializer=locremove.toDict())
			removeloc['removed'] = True
			self.locationToDB(removeloc)
		locnames = self.locationNames()
		self.panel.locationsEvent(locnames)

	def locationNames(self):
		names = [loc['name'] for loc in self.stagelocations]
		return names

	def getLocation(self, name):
		for i in self.stagelocations:
			if i['name'] == name:
				return i
		return None

	def getLocationNames(self):
		return map(lambda l: l['name'], self.stagelocations)

	def getLocationsFromDB(self):
		'''
		get list of locations for this session from DB
		and use them to create self.stagelocations list
		'''
		### get location from database
		locdata = leginondata.StageLocationData(session=self.session)
		locations = self.research(datainstance=locdata)

		### only want most recent of each name
		### since they are ordered by timestamp, just take the
		### first one of each name
		mostrecent = []
		names = []
		for loc in locations:
			if loc['name'] not in names:
				names.append(loc['name'])
				if not loc['removed']:
					mostrecent.append(loc)
		self.stagelocations[:] = mostrecent

	def locationToDB(self, stagelocdata):
		'''
		stores a location in the DB under the current session name
		if no location is specified, store all self.stagelocations
		'''
		self.publish(stagelocdata, database=True, dbforce=True)

	def toScope(self, loc):
		'''
		loc is either index, location, or name
		'''
		errstr = 'Set instrument failed: %s'
		locdata = None
		if type(loc) is int:
			locdata = self.stagelocations[p]
		elif type(loc) in (str, unicode):
			for location in self.stagelocations:
				if loc == location['name']:
					locdata = location
					break
		elif isinstance(loc, leginondata.StageLocationData):
			locdata = loc
		else:
			self.logger.error(errstr % 'bad argument')
			return

		if locdata is None:
			self.logger.error(errstr % 'no such location')
			return

		name = locdata['name']
		self.logger.info('Moving to location %s' % (name,))

		## should switch to using AllEMData
		stagedict = {}
		stagedict['x'] = locdata['x']
		stagedict['y'] = locdata['y']
		if not locdata['xy only']:
			stagedict['z'] = locdata['z']
			stagedict['a'] = locdata['a']
		self._toScope(name,stagedict)

	def _toScope(self,name, stagedict):
		try:
			self.instrument.tem.StagePosition = stagedict
		except KeyError:
			self.logger.exception('instrument key %s not available' % (stagedict.keys(),))
		except:
			self.logger.exception('unable to set instrument')
		else:
			self.logger.info('Moved to location %s' % (name,))

	def uiGetPreset(self,presetname):
		self.logger.info('Get %s preset image shift and beam shift from scope' % (presetname,))
		self.setStatus('processing')
		params = {}
		params['image shift'] = self.instrument.tem.ImageShift
		params['beam shift'] = self.instrument.tem.BeamShift
		self.presetsclient.updatePreset(presetname,params)
		self.setStatus('idle')

	def uiSendPreset(self,presetname):
		self.logger.info('Send %s preset to scope' % (presetname,))
		self.setStatus('processing')
		self.presetsclient.toScope(presetname)
		self.setStatus('idle')
		self.panel.onSendPresetDone()

	def onResetXY(self):
		loc = {'x':0.0,'y':0.0}
		self._toScope('reset xy',loc)

	def onResetZ(self):
		loc = {'z':0.0}
		self._toScope('reset z',loc)

	def onResetAlpha(self):
		loc = {'a':0.0}
		self._toScope('reset alpha',loc)

if __name__ == '__main__':
	id = ('navigator',)
	n = Navigator(id, None)


text:
import numpy as np

from ...bayesian_optimization.acquisitions import (
    ExpectedImprovement,
    NegativeLowerConfidenceBound,
    ProbabilityOfImprovement,
)
from ...core import ParameterSpace
from ...core.loop import FixedIntervalUpdater, OuterLoop, SequentialPointCalculator
from ...core.loop.loop_state import create_loop_state
from ...core.optimization import AcquisitionOptimizerBase, GradientAcquisitionOptimizer
from ..models.bohamiann import Bohamiann
from ..models.random_forest import RandomForest
from .enums import AcquisitionType, ModelType

def create_bayesian_optimization_loop(
    x_init: np.ndarray,
    y_init: np.ndarray,
    parameter_space: ParameterSpace,
    acquisition_type: AcquisitionType,
    model_type: ModelType,
    cost_init: np.ndarray = None,
    model_kwargs: dict = None,
    acquisition_optimizer: AcquisitionOptimizerBase = None,
) -> OuterLoop:
    """
    Creates Bayesian optimization loop for Bayesian neural network or random forest models.

    :param x_init: 2d numpy array of shape (no. points x no. input features) of initial X data
    :param y_init: 2d numpy array of shape (no. points x no. targets) of initial Y data
    :param cost_init: 2d numpy array of shape (no. points x no. targets) of initial cost of each function evaluation
    :param parameter_space: parameter space describing input domain
    :param acquisition_type: an AcquisitionType enumeration object describing which acquisition function to use
    :param model_type: A ModelType enumeration object describing which model to use.
    :param model_kwargs: Key work arguments for the model constructor. See individual models for details.
    :param acquisition_optimizer: Optimizer selecting next evaluation points by maximizing acquisition.
                                  Gradient based optimizer is used if None. Defaults to None.
    :return: OuterLoop instance
    """

    if model_kwargs is None:
        model_kwargs = dict()

    # Create model
    if model_type is ModelType.RandomForest:
        model = RandomForest(x_init, y_init, **model_kwargs)
    elif model_type is ModelType.BayesianNeuralNetwork:
        model = Bohamiann(x_init, y_init, **model_kwargs)
    else:
        raise ValueError("Unrecognised model type: " + str(model_type))

    # Create acquisition
    if acquisition_type is AcquisitionType.EI:
        acquisition = ExpectedImprovement(model)
    elif acquisition_type is AcquisitionType.PI:
        acquisition = ProbabilityOfImprovement(model)
    elif acquisition_type is AcquisitionType.NLCB:
        acquisition = NegativeLowerConfidenceBound(model)
    else:
        raise ValueError("Unrecognised acquisition type: " + str(acquisition_type))

    if acquisition_optimizer is None:
        acquisition_optimizer = GradientAcquisitionOptimizer(parameter_space)
    candidate_point_calculator = SequentialPointCalculator(acquisition, acquisition_optimizer)

    if cost_init is None:
        loop_state = create_loop_state(x_init, y_init)
    else:
        loop_state = create_loop_state(x_init, y_init, cost=cost_init)

    model_updater = FixedIntervalUpdater(model, 1)
    return OuterLoop(candidate_point_calculator, model_updater, loop_state)

ext: py

max_stars_repo_name: EmuKit/Emukit

max_issues_repo_name: EmuKit/Emukit

max_forks_repo_name: EmuKit/Emukit

reason: import numpy


text:
"""Command line tools for generating and plotting Fermi surfaces."""
import os
import sys
import warnings
from collections import defaultdict

import click
from click import option

from ifermi.defaults import AZIMUTH, ELEVATION, SCALE, SYMPREC, VECTOR_SPACING

plot_type = click.Choice(["matplotlib", "plotly", "mayavi"], case_sensitive=False)
spin_type = click.Choice(["up", "down"], case_sensitive=False)
property_type = click.Choice(["velocity", "spin"], case_sensitive=False)

# define shared Fermi surface generation options
_generation_options = [
    option("-f", "--filename", help="vasprun.xml file to plot"),
    option(
        "-m",
        "--mu",
        default=0.0,
        help="offset from the Fermi level at which to calculate Fermi surface",
        show_default=True,
    ),
    option(
        "--wigner/--no-wigner",
        "wigner_seitz",
        default=True,
        help="use Wigner-Seitz cell rather than reciprocal lattice parallelepiped",
        show_default=True,
    ),
    option(
        "-i",
        "--interpolation-factor",
        default=8.0,
        help="interpolation factor for band structure",
        show_default=True,
    ),
    option("--property", type=property_type, help="projection type"),
    option(
        "--projection-axis",
        nargs=3,
        type=float,
        help="use dot product of property onto cartesian axis (e.g. 0 0 1)",
    ),
    option(
        "--decimate-factor",
        type=float,
        help="factor by which to decimate surfaces (i.e. 0.8 gives 20 % fewer faces)",
    ),
    option(
        "--smooth", is_flag=True, help="smooth the Fermi surface", show_default=True
    ),
]

def generation_options(func):
    for opt in reversed(_generation_options):
        func = opt(func)
    return func

@click.group(
    context_settings=dict(help_option_names=["-h", "--help"], max_content_width=120)
)
def cli():
    """IFermi is a tool for the generation, analysis and plotting of Fermi surfaces."""

    def _warning(message, *_, **__):
        click.echo(f"WARNING: {message}\n")

    warnings.showwarning = _warning

    warnings.filterwarnings("ignore", category=UserWarning, module="matplotlib")
    warnings.filterwarnings("ignore", category=UnicodeWarning, module="matplotlib")
    warnings.filterwarnings("ignore", category=UserWarning, module="pymatgen")

@cli.command()
@generation_options
@option(
    "--norm/--no-norm",
    default=True,
    help="average property norms (overridden by --projection-axis)",
    show_default=True,
)
@option(
    "--precision",
    default=4,
    help="number of decimal places in output",
    show_default=True,
)
def info(filename, **kwargs):
    """Calculate information about the Fermi surface."""
    import numpy as np
    from tabulate import tabulate

    fs, bs = _get_fermi_surface(
        filename=filename,
        interpolation_factor=kwargs["interpolation_factor"],
        properties=kwargs["property"],
        mu=kwargs["mu"],
        decimate_factor=None,
        smooth=False,
        wigner_seitz=kwargs["wigner_seitz"],
        calculate_dimensionality=True,
    )

    np.set_printoptions(precision=kwargs["precision"])

    def dp(value):
        if isinstance(value, (list, tuple, np.ndarray)):
            return str(np.array(value))
        else:
            return f"{value:.{kwargs['precision']}g}"

    norm = kwargs["norm"]
    axis = kwargs["projection_axis"] or None
    prop = kwargs["property"]

    if prop == "velocity":
        unit = " m/s"
        table_unit = " [m/s]"
    else:
        unit = ""
        table_unit = ""
        prop = ""

    click.echo("Fermi Surface Summary\n=====================\n")
    click.echo(f"  # surfaces: {fs.n_surfaces}")
    click.echo(f"  Area: {dp(fs.area)} Å⁻²")

    if fs.has_properties:
        avg_property = fs.average_properties(norm=norm, projection_axis=axis)
        click.echo(f"  Avg {prop}: {dp(avg_property)}{unit}")

    for spin in fs.spins:
        if bs.is_spin_polarized:
            title = f"Spin {spin.name.capitalize()} Isosurfaces"
        else:
            title = "Isosurfaces"
        click.echo(f"\n{title}\n{len(title) * '~'}\n")

        table = defaultdict(list)
        for isosurface in fs.isosurfaces[spin]:
            table["Band"].append(isosurface.band_idx + 1)
            table["Area [Å⁻²]"].append(isosurface.area)

            avg_property = None
            if isosurface.has_properties:
                avg_property = isosurface.average_properties(norm, axis)
            table[f"{prop.capitalize()} avg{table_unit}"].append(avg_property)

            table["Dimensionality"].append(isosurface.dimensionality)
            table["Orientation"].append(isosurface.orientation)

        # filter columns that have None's in all rows
        table = {k: v for k, v in table.items() if any([i is not None for i in v])}

        # format table
        table_str = tabulate(
            table,
            headers=table.keys(),
            numalign="right",
            stralign="center",
            floatfmt=f"#.{kwargs['precision']}g",
        )

        # indent table 2 spaces in
        table_str = "  " + table_str.replace("\n", "\n  ")
        click.echo(table_str)

@cli.command()
@generation_options
@option("-o", "--output", "output_filename", help="output filename")
@option(
    "-s",
    "--symprec",
    default=SYMPREC,
    help="symmetry precision in Å",
    show_default=True,
)
@option(
    "-a",
    "--azimuth",
    default=AZIMUTH,
    help="viewpoint azimuth angle in °",
    show_default=True,
)
@option(
    "-e",
    "--elevation",
    default=ELEVATION,
    help="viewpoint elevation angle in °",
    show_default=True,
)
@option(
    "-t",
    "--type",
    "plot_type",
    default="plotly",
    type=plot_type,
    help="plotting type",
    show_default=True,
)
@option(
    "--color-property/--no-color-property",
    default=True,
    help="color Fermi surface properties",
    show_default=True,
)
@option("--property-colormap", help="matplotlib colormap name for properties")
@option(
    "--vector-property/--no-vector-property",
    help="show vector properties as arrows",
    show_default=True,
)
@option("--vector-colormap", help="matplotlib colormap name for vectors")
@option(
    "--vector-spacing",
    default=VECTOR_SPACING,
    help="spacing between property arrows",
    show_default=True,
)
@option("--cmin", type=float, help="minimum intensity on property colorbar")
@option("--cmax", type=float, help="maximum intensity on property colorbar")
@option("--vnorm", type=float, help="value by which to normalise vector lengths")
@option(
    "--scale-linewidth",
    is_flag=True,
    help="scale Fermi slice thickness by property",
    show_default=True,
)
@option(
    "--hide-surface", is_flag=True, help="hide the Fermi surface", show_default=True
)
@option(
    "--plot-index",
    multiple=True,
    default=None,
    help="plot specific bands indices (1 based)",
)
@option(
    "--hide-labels",
    is_flag=True,
    help="hide the high-symmetry k-point labels",
    show_default=True,
)
@option(
    "--hide-cell", is_flag=True, help="hide reciprocal cell boundary", show_default=True
)
@option("--spin", type=spin_type, help="select spin channel", show_default=True)
@option(
    "--slice",
    nargs=4,
    type=float,
    help="slice through the Brillouin zone (format: j k l dist)",
)
@option("--scale", default=SCALE, help="scale for image resolution", show_default=True)
def plot(filename, **kwargs):
    """Plot a Fermi surface from a vasprun.xml file."""
    from pymatgen.electronic_structure.core import Spin

    from ifermi.plot import FermiSlicePlotter, FermiSurfacePlotter, save_plot, show_plot

    try:
        import mayavi.mlab as mlab
    except ImportError:
        mlab = False

    if kwargs["property_colormap"]:
        kwargs["color_property"] = kwargs["property_colormap"]

    if kwargs["vector_colormap"]:
        kwargs["vector_property"] = kwargs["vector_colormap"]

    output_filename = kwargs["output_filename"]
    if mlab and kwargs["plot_type"] == "mayavi" and output_filename is not None:
        # handle mlab non interactive plots
        mlab.options.offscreen = True

    fs, _ = _get_fermi_surface(
        filename=filename,
        interpolation_factor=kwargs["interpolation_factor"],
        properties=kwargs["property"],
        mu=kwargs["mu"],
        decimate_factor=kwargs["decimate_factor"],
        smooth=kwargs["smooth"],
        wigner_seitz=kwargs["wigner_seitz"],
        calculate_dimensionality=False,
    )

    spin = {"up": Spin.up, "down": Spin.down, None: None}[kwargs["spin"]]
    projection_axis = kwargs["projection_axis"] or None

    if len(kwargs["plot_index"]) > 0:
        # user specified bands start at 1 but python api expects 0-based indices
        kwargs["plot_index"] = [int(i) - 1 for i in kwargs["plot_index"]]
    else:
        # click automatically converts plot_index to list, even though the default
        # should be None
        kwargs["plot_index"] = None

    if kwargs["slice"]:
        plane_normal = kwargs["slice"][:3]
        distance = kwargs["slice"][3]

        fermi_slice = fs.get_fermi_slice(plane_normal, distance)
        plotter = FermiSlicePlotter(fermi_slice, symprec=kwargs["symprec"])
        fig = plotter.get_plot(
            spin=spin,
            color_properties=kwargs["color_property"],
            vector_properties=kwargs["vector_property"],
            projection_axis=projection_axis,
            vector_spacing=kwargs["vector_spacing"],
            cmin=kwargs["cmin"],
            cmax=kwargs["cmax"],
            vnorm=kwargs["vnorm"],
            hide_slice=kwargs["hide_surface"],
            hide_labels=kwargs["hide_labels"],
            hide_cell=kwargs["hide_cell"],
            scale_linewidth=kwargs["scale_linewidth"],
            plot_index=kwargs["plot_index"],
        )
    else:
        plotter = FermiSurfacePlotter(fs, symprec=kwargs["symprec"])
        fig = plotter.get_plot(
            plot_type=kwargs["plot_type"],
            spin=spin,
            azimuth=kwargs["azimuth"],
            elevation=kwargs["elevation"],
            color_properties=kwargs["color_property"],
            vector_properties=kwargs["vector_property"],
            projection_axis=projection_axis,
            vector_spacing=kwargs["vector_spacing"],
            cmin=kwargs["cmin"],
            cmax=kwargs["cmax"],
            vnorm=kwargs["vnorm"],
            hide_surface=kwargs["hide_surface"],
            hide_labels=kwargs["hide_labels"],
            hide_cell=kwargs["hide_cell"],
            plot_index=kwargs["plot_index"],
        )

    if output_filename is None:
        show_plot(fig)
    else:
        click.echo(f"Saving plot to {output_filename}")
        save_plot(fig, output_filename, scale=kwargs["scale"])

def find_vasprun_file():
    """Search for vasprun files from the current directory.

    Will look for vasprun.xml or vasprun.xml.gz files.
    """
    for file in ["vasprun.xml", "vasprun.xml.gz"]:
        if os.path.exists(file):
            return file

    click.echo("ERROR: No vasprun.xml found in current directory")
    sys.exit()

def _get_fermi_surface(
    filename,
    interpolation_factor,
    properties,
    mu,
    decimate_factor,
    smooth,
    wigner_seitz,
    calculate_dimensionality,
):
    """Common helper method to get Fermi surface"""
    import numpy as np
    from pymatgen.electronic_structure.core import Spin
    from pymatgen.io.vasp.outputs import Vasprun

    from ifermi.interpolate import FourierInterpolator
    from ifermi.kpoints import kpoints_from_bandstructure
    from ifermi.surface import FermiSurface

    if not filename:
        filename = find_vasprun_file()

    parse_projections = properties == "spin"
    vr = Vasprun(filename, parse_projected_eigen=parse_projections)
    bs = vr.get_band_structure()

    interpolator = FourierInterpolator(bs)
    interp_bs, velocities = interpolator.interpolate_bands(
        interpolation_factor, return_velocities=True
    )

    property_data = None
    property_kpoints = None
    if properties == "velocity":
        property_data = velocities
        property_kpoints = kpoints_from_bandstructure(interp_bs)
    elif properties == "spin":
        if vr.projected_magnetisation is not None:
            # transpose so shape is (nbands, nkpoints, natoms, norbitals, 3)
            property_data = vr.projected_magnetisation.transpose(1, 0, 2, 3, 4)

            # sum across all atoms and orbitals
            property_data = property_data.sum(axis=(2, 3))
            property_data /= np.linalg.norm(property_data, axis=-1)[..., None]

            property_data = {Spin.up: property_data}
            property_kpoints = kpoints_from_bandstructure(bs)
        else:
            click.echo(
                "ERROR: Band structure does not include spin properties.\n"
                "Ensure calculation was run with LSORBIT or LNONCOLLINEAR = True "
                "and LSORBIT = 11."
            )
            sys.exit()

    return (
        FermiSurface.from_band_structure(
            interp_bs,
            mu=mu,
            wigner_seitz=wigner_seitz,
            decimate_factor=decimate_factor,
            smooth=smooth,
            property_data=property_data,
            property_kpoints=property_kpoints,
            calculate_dimensionality=calculate_dimensionality,
        ),
        bs,
    )

max_stars_repo_name: whzhangg/IFermi

max_issues_repo_name: skphy/IFermi

max_forks_repo_name: skphy/IFermi


text:
'''
Pure python implementation of a connect 4 terminal game object.
Optimizations applied allow computation of one move and one check in approx. 100us.
---Still kinda slow... see connect4tf.py for a (hopefully) faster implementation.---nvm, this is all i got
'''
import numpy as np

class Connect4Board(object):
    def __init__(self, board_shape=(6, 7), winVecs=None):
        if winVecs is not None:
            self.wins=winVecs
        self.grid = np.zeros(board_shape, dtype=np.int8)
        self.height = np.zeros(board_shape[1], dtype=np.int8)
        self.player = 1
    def move(self, prob):# Add a piece to the board
        prob = prob * np.less(self.height, 6)
        if np.sum(prob) > 0:
            slot = np.argmax(prob)
            self.grid[self.height[slot], slot] = self.player
            self.height[slot] += 1
            #self.player *= -1# this is a bit slower than if statements...
            if self.player < 0:# swap players
                self.player = 1
            else:
                self.player = -1
            return None# No tie
        return 0# Game is a tie
    def check(self):# Check if anyone has won
        checked = np.dot(self.grid.reshape(42,), self.wins)
        if np.sum(checked > 3):
            return 1
        if np.sum(checked < -3):
            return -1
        return None# No winner
    def __str__(self):# String for print
        return str(self.grid[::-1])

class BoardExplorer(Connect4Board):
    def __init__(self,  board_shape=(6, 7), toWin=4):
        super().__init__(board_shape)
        self.toWin = toWin
        self.wins = []
    def findHorizWins(self):
        for i in range(self.grid.shape[1]-self.toWin+1):
            horiz = np.zeros(self.grid.shape[1], dtype=np.int8)
            for k in range(self.toWin):
                horiz[i+k] = 1
            for j in range(self.grid.shape[0]):
                self.grid[j] = horiz
                self.wins.append(self.grid)
                # self.wins.append(np.where(self.grid > 0))
                # print(super().__str__())
                super().__init__(self.grid.shape) #reset grid
        # print("winpatterns:", len(self.wins))
    def findVerticalWns(self):
        for i in range(self.grid.shape[0]-self.toWin+1):
            vert = np.zeros((self.grid.shape[0], 1), dtype=np.int8)
            for k in range(self.toWin):
                vert[i+k][0] = 1
            for j in range(self.grid.shape[1]):
                self.grid[:,j] = vert[:,0]
                self.wins.append(self.grid)
                # self.wins.append(np.where(self.grid > 0))
                # print(super().__str__())
                super().__init__(self.grid.shape) #reset grid
        # print("winpatterns:", len(self.wins))
    def findDiagWins(self):
        # first half:
        row = 0
        while(row < self.grid.shape[0]):
            col = 0
            rowTmp = row
            diag = []
            while rowTmp >= 0: # get the diagonal
                diag.append((rowTmp, col))
                rowTmp -= 1
                col += 1
            row += 1
            if len(diag) >= self.toWin:# if diag large enough,
                # step through all positions of a winning sequence on the grid:
                for i in range(len(diag)-self.toWin+1):
                    # Choose points:
                    winDiag = [diag[i+k] for k in range(self.toWin)]
                    # convert points to np.array indices
                    winIdx = (np.array([i for i, j in winDiag]), np.array([j for i, j in winDiag]))
                    self.grid[winIdx] = 1# winning diag = to 1
                    self.wins.append(self.grid)
                    # self.wins.append(np.where(self.grid > 0))#save points
                    # print(super().__str__())
                    self.grid = np.flip(self.grid, axis=0)# flip vertically
                    self.wins.append(self.grid)
                    # self.wins.append(np.where(self.grid > 0))#save points
                    # print(super().__str__())
                    self.grid = np.flip(self.grid, axis=1)# flip horizontally
                    self.wins.append(self.grid)
                    # self.wins.append(np.where(self.grid > 0))#save points
                    # print(super().__str__())
                    self.grid = np.flip(self.grid, axis=0)# flip vertically
                    self.wins.append(self.grid)
                    # self.wins.append(np.where(self.grid > 0))#save points
                    # print(super().__str__())
                    super().__init__(self.grid.shape)# reset grid
        print("winpatterns:", len(self.wins))
    def getWinPatterns(self):
        self.wins = []
        self.findHorizWins()
        self.findVerticalWns()
        self.findDiagWins()
        filters = np.swapaxes(np.array(self.wins, dtype=np.int8).T, 0, 1) #(6x7x69)
        return filters.reshape(42, 69) # single vector for comparison

if __name__ == "__main__":
    explore = BoardExplorer()
    # explore.findHorizWins()
    # explore.findVerticalWns()
    # explore.findDiagWins()
    winVecs = explore.getWinPatterns()
    # np.save("winVecs", winVecs)
    board = Connect4Board(winVecs=winVecs)
    for i in range(5):
        board.move(i+1)
        board.move(i)
        board.move(i+1)
        board.move(i)
        board.move(i+1)
        board.move(i)
    print(board)
    print(board.check())

    winVecs1 = explore.getWinPatterns()

    board = Connect4Board(winVecs=winVecs1)
    import timeit
    print(timeit.timeit("board = Connect4Board(winVecs=winVecs)", setup="from __main__ import Connect4Board, winVecs", number=1000)/1000)
    print(timeit.timeit("board.move(1);board = Connect4Board(winVecs=winVecs)", setup="from __main__ import board, Connect4Board, winVecs", number=10000)/10000)
    print(timeit.timeit("board.check()", setup="from __main__ import board", number=10000)/10000)

    # x = board.grid.reshape((42,))
    # y = winFilters
    # print(np.dot(x, y))
    # exit()

    # print(winFilters)
    # print(winFilters.shape)
    # x = np.matmul(board.grid, winFilters)
    # print(x[:,:,:])
    # x = np.sum(x, axis=0)
    # x = np.sum(x, axis=0)
    # # print(winFilters[:,:,np.where(x == 16)[0][0]])
    # print(x.shape)
    # print(x)
    # exit()

    # print("wins = [")
    # for i in wins:
    #     print("(np.array([{},{},{},{}], dtype=np.int32), np.array([{},{},{},{}], dtype=np.int32)),".format(
    #     i[0][0], i[0][1], i[0][2], i[0][3], i[1][0], i[1][1], i[1][2], i[1][3]))
    # print("]")
    # exit()

    # import timeit
    # x = np.ndarray((6, 7), dtype=np.int8)
    # y = np.ndarray((6, 7, 69), dtype=np.int8)
    # print(timeit.timeit("np.sum(np.dot(x, y))", setup="from __main__ import x, y, np", number=10000)/10000)
    # print(timeit.timeit("np.sum(np.dot(board, winFilters))", setup="from __main__ import board, winFilters, np", number=1000)/1000)
    # print(timeit.timeit("board.check()", setup="from __main__ import board", number=1000)/1000)

    # print(timeit.timeit("board = Connect4Board()", setup="from __main__ import Connect4Board", number=1000)/1000)
    # import time
    # N = 1000
    # start = time.time()
    # for i in range(N):
    #     board = Connect4Board()
    #     for k in range(6):
    #         for j in range(7):
    #             board.move(j)
    #             [np.sum(board.grid[i]) for i in wins]
    #
    # print(1000/(time.time()-start))


text:
#! /usr/bin/env python
# coding=utf-8
#================================================================
#   Copyright (C) 2019 * Ltd. All rights reserved.
#
#   Editor      : VIM
#   File name   : fuck.py
#   Author      : YunYang1994
#   Created date: 2019-01-23 10:21:50
#   Description :
#
#================================================================

import numpy as np
import tensorflow as tf
from PIL import Image
from core import utils
import os
os.environ["CUDA_DEVICE_ORDER"]="PCI_BUS_ID"
os.environ["CUDA_VISIBLE_DEVICES"] = str(0)

IMAGE_H, IMAGE_W = 416, 416
classes = utils.read_coco_names('./sixray.names')     # todo
num_classes = len(classes)
image_path = "/media/kent/DISK2/SBRI_Project/dataset_sixray/SIXray-master/ks3util-1.1.1-upload/JPEGImage/P08852.jpg"  # 181,
img = Image.open(image_path)
img_resized = np.array(img.resize(size=(IMAGE_W, IMAGE_H)), dtype=np.float32)
img_resized = img_resized / 255.
cpu_nms_graph = tf.Graph()

input_tensor, output_tensors = utils.read_pb_return_tensors(cpu_nms_graph, "./checkpoint/yolov3_cpu_nms.pb",
                                           ["Placeholder:0", "concat_9:0", "mul_6:0"])
with tf.Session(graph=cpu_nms_graph) as sess:
    boxes, scores = sess.run(output_tensors, feed_dict={input_tensor: np.expand_dims(img_resized, axis=0)})
    boxes, scores, labels = utils.cpu_nms(boxes, scores, num_classes, score_thresh=0.3, iou_thresh=0.5)
    image = utils.draw_boxes(img, boxes, scores, labels, classes, [IMAGE_H, IMAGE_W], show=True)


text:
import numpy as np
import sympy as sp
from sympy.solvers.solveset import linsolve
from pprint import pprint

def get_cycle(basis_indexes, start, finish):
    used = [False for i in range(len(basis_indexes))]

    def dfs(basis_indexes, start, finish):
        straight = []
        reverse = []

        if start == finish:
            return straight, reverse

        for i, (from_, to_) in enumerate(basis_indexes):
            if used[i]:
                continue

            if from_ == start:
                used[i] = True
                res = dfs(basis_indexes, to_, finish)
                if res:
                    straight.append((from_, to_))
                    straight += res[0]
                    reverse += res[1]
                    return straight, reverse

            if to_ == start:
                used[i] = True
                res = dfs(basis_indexes, from_, finish)
                if res:
                    reverse.append((from_, to_))
                    straight += res[0]
                    reverse += res[1]
                    return straight, reverse

    straight, reverse = dfs(basis_indexes, start, finish)
    straight.append((finish, start))

    return straight, reverse

def get_c(costs, i, j):
    for to, cost in costs[i]:
        if to == j:
            return cost

def get_u(costs, basis):
    n = len(costs)
    u = np.zeros(n)
    used = [False for _ in range(n)]
    used[0] = True

    while not all(used):
        for i, j in basis:
            if used[i]:
                u[j] = u[i] - get_c(costs, i, j)
                used[j] = True
            elif used[j]:
                u[i] = u[j] + get_c(costs, i, j)
                used[i] = True

    return u

def min_cost_edge(edges):
    min_x = 1e9
    min_pair = None

    for from_, to_ in edges:
        cost = get_c(x, from_, to_)
        if cost < min_x:
            min_x = cost
            min_pair = (from_, to_)

    return min_x, min_pair

def get_not_basis(graph, basis):
    not_basis_edges = []
    for src, dests in enumerate(graph):
        for dest in dests:
            if (src, dest) not in basis:
                not_basis_edges.append((src, dest))
    return not_basis_edges

def get_delta(costs, u, not_basis_edges):
    delta = []
    for i, j in not_basis_edges:
        delta.append(u[i] - u[j] - get_c(costs, i, j))
    return delta

def get_bad_edge(costs, u, not_basis_edges):
    for i, j in not_basis_edges:
        if u[i] - u[j] - get_c(costs, i, j) > 0:
            return (i, j)

    return None

def run(costs, basis_indexes, x):
    graph = []
    for d in costs:
        g = []
        for e in d:
            g.append(e[0])
        graph.append(g)

    print('X array:', x)
    print('Graph:', graph)
    print('b', basis_indexes)
    print()

    count = 0
    while True:
        count += 1
        print('\nIter', count)

        u = get_u(costs, basis_indexes)
        print('u', u)

        not_basis_edges = get_not_basis(graph, basis_indexes)
        print('Not basis:', not_basis_edges)

        delta = get_delta(costs, u, not_basis_edges)
        print('Estimate vector:', delta)

        bad = get_bad_edge(costs, u, not_basis_edges)
        if not bad:
            return x, basis_indexes

        print('Bad is:', bad)

        straight, reverse = get_cycle(basis_indexes, bad[1], bad[0])
        print('Straight', straight)
        print('Reverse', reverse)

        min_x, min_pair = min_cost_edge(reverse)

        for from_, to_ in straight:
            for i in range(len(x[from_])):
                if x[from_][i][0] == to_:
                    x[from_][i] = (x[from_][i][0],
                                  x[from_][i][1] + min_x)

        for r in reverse:
            for i in range(len(x[r[0]])):
                if x[r[0]][i][0] == r[1]:
                    x[r[0]][i] = (x[r[0]][i][0],
                                  x[r[0]][i][1] - min_x)

        basis_indexes.remove(min_pair)
        basis_indexes.append(bad)

def get_data():
    data = [
        [(1, 1)],
        [(5, 3)],
        [(1, 3), (3, 5)],
        [],
        [(2, 4), (3, 1)],
        [(0, -2), (2, 3), (4, 4)]
    ]
    basis = [
        (0, 1),
        (2, 1),
        (5, 2),
        (2, 3),
        (4, 3)
    ]
    x = [
        [(1, 1)],
        [(5, 0)],
        [(1, 3), (3, 1)],
        [],
        [(2, 0), (3, 5)],
        [(0, 0), (2, 9), (4, 0)]
    ]
    return data, basis, x

def get_sample():
    data = [
        [(1, 9), (7, 5)],
        [(2, 1), (5, 3), (6, 5)],
        [(8, -2)],
        [(2, -3)],
        [(3, 6)],
        [(4, 8)],
        [(2, -1), (3, 4), (4, 7), (8, 1)],
        [(6, 2), (8, 2)],
        [(5, 6)]
    ]
    basis = [
        (0, 1),
        (0, 7),
        (1, 6),
        (1, 2),
        (6, 4),
        (4, 3),
        (5, 4),
        (8, 5)
    ]
    x = [
        [(1, 2), (7, 7)],
        [(2, 4), (5, 0), (6, 3)],
        [(8, 0)],
        [(2, 0)],
        [(3, 3)],
        [(4, 4)],
        [(2, 0), (3, 0), (4, 5), (8, 0)],
        [(6, 0), (8, 0)],
        [(5, 2)]
    ]
    return data, basis, x

data, basis_indexes, x = get_data()
answer = run(data, basis_indexes, x)
if answer is not None:
    print('Answer. X: {}'.format(answer[0]))
    x = answer[0]
    d = data
    s = 0
    for i in range(len(x)):
        for j in range(len(x[i])):
            s += x[i][j][1] * get_c(d, i, x[i][j][0])
    print('S:', s)

max_stars_repo_name: andrew-kulikov/saiio

max_issues_repo_name: andrew-kulikov/saiio

max_forks_repo_name: andrew-kulikov/saiio

reason: import numpy,import sympy,from sympy


text:
#!/usr/bin/env python
from __future__ import print_function

import torch
from torch import nn, optim
import torch.nn.init as init

import cv2
import sys
sys.path.append("game/")

import wrapped_flappy_bird as game
import random
import numpy as np
from collections import deque

GAME = 'bird' # the name of the game being played for log files
ACTIONS = 2 # number of valid actions
GAMMA = 0.99 # decay rate of past observations
OBSERVE = 1370001. # timesteps to observe before training
EXPLORE = 2000000. # frames over which to anneal epsilon
FINAL_EPSILON = 0.001 # final value of epsilon
INITIAL_EPSILON = 0.1 # starting value of epsilon
REPLAY_MEMORY = 50000 # number of previous transitions to remember
BATCH = 64 # size of minibatch
FRAME_PER_ACTION = 1

torch.set_default_tensor_type('torch.cuda.FloatTensor')
device = torch.device('cuda')
# batch_size = 64
# image_size = 80

# #Create transform
# transform = transforms.Compose(
#     [transforms.Resize(image_size),
#     transforms.ToTensor(),
#     transforms.Normalize((0.5,0.5,0.5),(0.5,0.5,0.5)),]
# )

# #Load dataset
# dataset = datasets.CIFAR10(root = 'F:/t/python/pytorch-practice/cifar-10', download=False,transform=transform)
# dataloader = torch.utils.data.DataLoader(dataset, batch_size= batch_size, shuffle=True, num_workers=2)
def weight_init(m):
    classname = m.__class__.__name__
    if classname.find('Conv') != -1:
        m.weight.data.normal_(0.0, 1.0)
        m.bias.data.fill_(0.1)
    elif classname.find('Linear') != -1:
        m.weight.data.normal_(0.0,1.0)
        m.bias.data.fill_(0.1)

class Net(nn.Module):
    def __init__(self):
        super(Net, self).__init__()
        self.conv = nn.Sequential(
            nn.Conv2d(4, 32, 8, bias=True, stride=4,padding=2),
            #nn.ReLU(inplace=True),
            nn.MaxPool2d(2,stride=2),
            nn.Conv2d(32,64,4, bias= True,stride=2,padding=1),
            #nn.ReLU(inplace=True),
            nn.MaxPool2d(2,stride=2,padding=1),
            nn.Conv2d(64,64,3,bias=True,stride=1, padding=1),
            #nn.ReLU(inplace=True),
            nn.MaxPool2d(2,stride=1)
        )
        self.fc = nn.Sequential(
            nn.ReLU(),
            nn.Linear(256, 2,bias=True)
            #nn.Linear(512,ACTIONS,bias=True),
        )
        for layer in self.conv.children() :
            if isinstance(layer, nn.Conv2d):
                init.normal_(layer.weight, mean = 0, std=0.01)
                init.constant_(layer.bias, 0.1)

        for layer in self.fc.children():
            if isinstance(layer, nn.Linear):
                init.normal_(layer.weight,std=0.01)
                init.constant_(layer.bias, 0.1)

        #torch.autograd.set_grad_enabled(True)

    def forward(self, inputs):
        return self.fc(self.conv(inputs).view(-1,256))

# define the cost function
def cost(a, readout, y):
    readout_action = (readout* a).sum(dim = 1)
    t = y.detach()- readout_action
    cost = torch.mean(t*t)
    return cost

def trainNetwork(net):
    loss_func = nn.MSELoss()
    optimizer = optim.Adam(net.parameters(), lr = 1e-6)

    # open up a game state to communicate with emulator
    game_state = game.GameState()

    # store the previous observations in replay memory
    D = deque()

    # printing
    # a_file = open("logs_" + GAME + "/readout.txt", 'w')
    # h_file = open("logs_" + GAME + "/hidden.txt", 'w')

    # get the first state by doing nothing and preprocess the image to 80x80x4
    do_nothing = np.zeros(ACTIONS)
    do_nothing[0] = 1
    x_t, r_0, terminal = game_state.frame_step(do_nothing)
    current_frame = cv2.cvtColor(cv2.resize(x_t, (80, 80)), cv2.COLOR_BGR2GRAY)
    ret,current_frame = cv2.threshold(current_frame,1,255,cv2.THRESH_BINARY)
    current_frame.astype(np.float32)
    frame_data = np.stack((current_frame, current_frame, current_frame, current_frame), axis=0)
    input_img = torch.from_numpy(frame_data).float()

    # start training
    epsilon = 0.0001
    t = 0
    while "flappy bird" != "angry bird":
        # choose an action epsilon greedily
        Q_value1 = net(input_img.cuda().unsqueeze(0))
        a_t = np.zeros([ACTIONS])

        if t % FRAME_PER_ACTION == 0 :
            if random.random() <= epsilon:
                print("----------Random Action----------")
                action_index = random.randrange(ACTIONS)
                a_t[action_index] = 1
            else:
                action_index = np.argmax(Q_value1.cpu().detach().numpy())
                a_t[action_index] = 1
        else:
            a_t[0] = 1 # do nothing

        # scale down epsilon
        if epsilon > FINAL_EPSILON and t > OBSERVE:
            epsilon -= (INITIAL_EPSILON - FINAL_EPSILON) / EXPLORE

        # run the selected action and observe next state and reward
        x_t1_colored, r_t, terminal = game_state.frame_step(a_t)
        frame = cv2.cvtColor(cv2.resize(x_t1_colored, (80, 80)), cv2.COLOR_BGR2GRAY)
        ret,frame = cv2.threshold(frame,1,255,cv2.THRESH_BINARY)
        frame.astype(np.float32)
        frame = torch.Tensor(frame).cpu()
        new_input = torch.stack([frame, input_img[0,:,:],input_img[1,:,:],input_img[2,:,:]])
        #new_frame_data = np.stack((frame, current_frame, current_frame, current_frame), axis=0)
        #new_input = torch.from_numpy(new_frame_data).float()

        # store the transition in D
        D.append(
            (input_img,
            torch.Tensor(a_t, device = device),
            torch.Tensor([r_t], device = device),
            new_input,
            terminal)
        )
        if len(D) > REPLAY_MEMORY:
            D.popleft()

        # only train if done observing
        if t > OBSERVE:
            # sample a minibatch to train on
            minibatch = random.sample(D, BATCH)

            # get the batch variables
            orginal_input_batch =  torch.stack([d[0] for d in minibatch]).cuda()#s_t
            a_batch = torch.stack([d[1] for d in minibatch]).cpu()#a_t
            reward_batch = torch.stack([d[2] for d in minibatch]).cuda()#r_t
            current_input_batch = torch.stack([d[3] for d in minibatch]).cuda()#s_t1

            action_batch = a_batch.numpy()
            index = np.argmax(action_batch, axis=1)
            index=np.reshape(index,[BATCH,1])
            index = torch.Tensor(index).cuda()
            y_batch = []
            Q_value_batch = net(current_input_batch)
            # for inputs in current_input_batch:
            #     readout_j1_batch.append(net(inputs.cuda()))#GPU

            for i in range(0, len(minibatch)):
                terminal = minibatch[i][4]
                # if terminal, only equals reward
                if terminal:
                    y_batch.append(reward_batch[i])
                else:
                    y_batch.append(reward_batch[i] + GAMMA * torch.max(Q_value_batch[i].detach()))

            y_batch = torch.stack(y_batch).cuda()
            y_predict = net(orginal_input_batch).gather(1,index.long())
            #y_predict.requires_grad_()
            # perform gradient step

            # for i in range(0, len(y_batch)):
            loss = loss_func(y_predict,y_batch)
            #print( 'loss: ',loss)
            optimizer.zero_grad()
            loss.backward()
            optimizer.step()

        # update the old values
        input_img = new_input
        t += 1

        # save progress every 10000 iterations
        if t % 10000 == 0:
            torch.save(net,'dqn-new.pkl')
            #saver.save(sess, 'saved_networks/' + GAME + '-dqn', global_step = t)

        # print info
        state = ''
        if t <= OBSERVE:
            state = "observe"
        elif t > OBSERVE and t <= OBSERVE + EXPLORE:
            state = "explore"
        else:
            state = "train"

        print("TIMESTEP", t, "/ STATE", state, \
            "/ EPSILON", epsilon, "/ ACTION", action_index, "/ REWARD", r_t, \
            "/ Q_MAX %e" % np.max(Q_value1.cpu().detach().numpy()))

if __name__ == "__main__":
    #net = Net().cuda()
    net = torch.load('dqn-2.pkl')
    #torch.autograd.set_grad_enabled(True)
    #net = net.apply(weight_init)
    trainNetwork(net)


text:
from glob import glob
import pandas as pd
import sys
import ntpath
import numpy as np
from PIL import Image, ImageDraw, ImageFont
from nltk.tokenize import RegexpTokenizer
from skimage.transform import resize
import os
import errno

checkpoint = '2018_10_23_14_55_12'
gen_dir = '../data/coco/gen_masks_%s/*'%(checkpoint)
gt_dir = '../data/coco/masks/'
img_dir = '../data/coco/images/'
txt_dir = '../data/coco/text/'
output_dir = '../vis_bboxes_%s/'%(checkpoint)
CAPS_PER_IMG = 5
FONT_MAX = 40
FONT_REAL = 30
MAX_WORD_NUM = 20
FNT = ImageFont.truetype('../data/coco/share/Pillow/Tests/fonts/FreeMono.ttf', FONT_REAL)
STD_IMG_SIZE = 256
VIS_SIZE = STD_IMG_SIZE
OFFSET = 2
SHOW_LIMIT = 500

def path_leaf(path):
	return ntpath.basename(path)

def load_captions(cap_path):
    all_captions = []
    with open(cap_path, "r") as f:
        captions = f.read().decode('utf8').split('\n')
        cnt = 0
        for cap in captions:
            if len(cap) == 0:
                continue
            cap = cap.replace("\ufffd\ufffd", " ")
            # picks out sequences of alphanumeric characters as tokens
            # and drops everything else
            tokenizer = RegexpTokenizer(r'\w+')
            tokens = tokenizer.tokenize(cap.lower())
            # print('tokens', tokens)
            if len(tokens) == 0:
                print('cap', cap)
                continue

            tokens_new = []
            for t in tokens:
                t = t.encode('ascii', 'ignore').decode('ascii')
                if len(t) > 0:
                    tokens_new.append(t)
            sentence = ' '.join(tokens_new)
            all_captions.append(sentence)
            cnt += 1
            if cnt == CAPS_PER_IMG:
                break
        if cnt < CAPS_PER_IMG:
            print('ERROR: the captions for %s less than %d'
                  % (cap_path, cnt))
    return all_captions

def draw_plate(bboxes):
	bbox_plate = Image.fromarray((np.ones((VIS_SIZE, VIS_SIZE, 3))*255).astype(np.uint8))
	if bboxes is None:
		return bbox_plate

	d = ImageDraw.Draw(bbox_plate)
	for i in xrange(bboxes.shape[0]):
		left, top, width, height, label = bboxes[i, :5]
		label = int(label)
		color = (210-label*2,label*3,50+label*2)
		d.rectangle([left, top, left+width-1, top+height-1], outline=color)
		d.text([left+5, top+5], str(label), fill=color)
	del d
	return bbox_plate

def mkdir_p(path):
    try:
        os.makedirs(path)
    except OSError as exc:  # Python >2.5
        if exc.errno == errno.EEXIST and os.path.isdir(path):
            pass
        else:
            raise

def is_non_zero_file(fpath):
    return True if os.path.isfile(fpath) and os.path.getsize(fpath) > 0 else False

mkdir_p(output_dir)

gen_paths = glob(gen_dir)
keys = [path_leaf(gen_path) for gen_path in gen_paths]

count = 0
for key in keys:
	if count >= SHOW_LIMIT:
		break

	# 1. load image
	img_path = '%s%s.jpg'%(img_dir, key)
	img = np.array(Image.open(img_path))
	img_height, img_width = img.shape[0], img.shape[1]
	height_scale = STD_IMG_SIZE/float(img_height)
	width_scale = STD_IMG_SIZE/float(img_width)
	img = resize(img, [STD_IMG_SIZE, STD_IMG_SIZE])
	img = Image.fromarray((img*255).astype(np.uint8))

	# 2. load captions
	cap_path = '%s%s.txt'%(txt_dir, key)
	captions = load_captions(cap_path)

	# 3. load gt bboxes
	gt_bbox_path = '%s%s/boxes.txt'%(gt_dir, key)
	if is_non_zero_file(gt_bbox_path):
		gt_boxes = pd.read_csv(gt_bbox_path, header=None).astype(int)
		gt_boxes = np.array(gt_boxes)
		gt_boxes[:,[0,2]] = gt_boxes[:,[0,2]]*width_scale
		gt_boxes[:,[1,3]] = gt_boxes[:,[1,3]]* height_scale
	else:
		gt_boxes = None

	# 4. load gen bboxes
	gen_bbox_paths = glob('%s%s/*'%(gen_dir, key))
	gen_bbox_paths_indices = [int(path_leaf(gen_bbox_path)) for gen_bbox_path in gen_bbox_paths]
	gen_bbox_paths_indices = np.argsort(gen_bbox_paths_indices)
	gen_bbox_paths = [gen_bbox_paths[index] for index in gen_bbox_paths_indices]
	gen_boxes_set = []
	for gen_bbox_path in gen_bbox_paths:
		sub_gen_bbox_path = '%s/boxes.txt'%(gen_bbox_path)
		if is_non_zero_file(sub_gen_bbox_path):
			gen_boxes = pd.read_csv(sub_gen_bbox_path, header=None).astype(int)
			gen_boxes = np.array(gen_boxes)
			gen_boxes_set.append(gen_boxes)
		else:
			gen_boxes_set.append(None)

	# 5. draw text
	text_convas = np.ones([VIS_SIZE + CAPS_PER_IMG * FONT_MAX,  7 * VIS_SIZE, 3], dtype=np.uint8)
	img_plate = Image.fromarray(text_convas)

	d = ImageDraw.Draw(img_plate)
	for i in xrange(len(captions)):
		d.text((2, i*FONT_REAL), captions[i], font=FNT, fill=(255, 255, 255, 255))

	del d

	# 6. draw images
	starting_y = CAPS_PER_IMG * FONT_MAX
	starting_x = VIS_SIZE+OFFSET

	img_plate.paste(img, (0, starting_y))

	# 7. draw gt boxes
	gt_bbox_plate = draw_plate(gt_boxes)
	img_plate.paste(gt_bbox_plate, (starting_x, starting_y))

	# 8. draw gen boxes
	for i in xrange(len(gen_boxes_set)):
		gen_boxes = gen_boxes_set[i]
		gen_bbox_plate = draw_plate(gen_boxes)
		img_plate.paste(gen_bbox_plate, (starting_x*(i+2), starting_y))

	img_plate.save('%s%s.jpg'%(output_dir, key))


text:
"""Utilities for norms and metrics."""
from typing import Any, Callable, Generic, Optional, Tuple, TypeVar

import multimethod
import numpy as np

from ..._utils import _pairwise_symmetric
from ...representation import FData, FDataGrid
from ...representation._typing import NDArrayFloat, Vector
from ._typing import Metric, MetricElementType, Norm, VectorType

T = TypeVar("T", bound=FData)

def _check_compatible(fdata1: T, fdata2: T) -> None:

    if isinstance(fdata1, FData) and isinstance(fdata2, FData):
        if (
            fdata2.dim_codomain != fdata1.dim_codomain
            or fdata2.dim_domain != fdata1.dim_domain
        ):
            raise ValueError("Objects should have the same dimensions")

        if not np.array_equal(fdata1.domain_range, fdata2.domain_range):
            raise ValueError("Domain ranges for both objects must be equal")

def _cast_to_grid(
    fdata1: FData,
    fdata2: FData,
    eval_points: Optional[NDArrayFloat] = None,
    _check: bool = True,
) -> Tuple[FDataGrid, FDataGrid]:
    """
    Convert fdata1 and fdata2 to FDatagrid.

    Checks if the fdatas passed as argument are unidimensional and compatible
    and converts them to FDatagrid to compute their distances.

    Args:
        fdata1: First functional object.
        fdata2: Second functional object.
        eval_points: Evaluation points.

    Returns:
        Tuple with two :obj:`FDataGrid` with the same grid points.
    """
    # Dont perform any check
    if not _check:
        return fdata1, fdata2

    _check_compatible(fdata1, fdata2)

    # Case new evaluation points specified
    if eval_points is not None:  # noqa: WPS223
        fdata1 = fdata1.to_grid(eval_points)
        fdata2 = fdata2.to_grid(eval_points)

    elif not isinstance(fdata1, FDataGrid) and isinstance(fdata2, FDataGrid):
        fdata1 = fdata1.to_grid(fdata2.grid_points[0])

    elif not isinstance(fdata2, FDataGrid) and isinstance(fdata1, FDataGrid):
        fdata2 = fdata2.to_grid(fdata1.grid_points[0])

    elif (
        not isinstance(fdata1, FDataGrid)
        and not isinstance(fdata2, FDataGrid)
    ):
        domain = fdata1.domain_range[0]
        grid_points = np.linspace(*domain)
        fdata1 = fdata1.to_grid(grid_points)
        fdata2 = fdata2.to_grid(grid_points)

    elif not np.array_equal(
        fdata1.grid_points,
        fdata2.grid_points,
    ):
        raise ValueError(
            "Grid points for both objects must be equal or"
            "a new list evaluation points must be specified",
        )

    return fdata1, fdata2

class NormInducedMetric(Metric[VectorType]):
    r"""
    Metric induced by a norm.

    Given a norm :math:`\| \cdot \|: X \rightarrow \mathbb{R}`,
    returns the metric :math:`d: X \times X \rightarrow \mathbb{R}` induced
    by the norm:

    .. math::
        d(f,g) = \|f - g\|

    Args:
        norm: Norm used to induce the metric.

    Examples:
        Computes the :math:`\mathbb{L}^2` distance between an object containing
        functional data corresponding to the function :math:`y(x) = x` defined
        over the interval [0, 1] and another one containing data of the
        function :math:`y(x) = x/2`.

        Firstly we create the functional data.

        >>> import skfda
        >>> import numpy as np
        >>> from skfda.misc.metrics import l2_norm, NormInducedMetric
        >>>
        >>> x = np.linspace(0, 1, 1001)
        >>> fd = skfda.FDataGrid([x], x)
        >>> fd2 = skfda.FDataGrid([x/2], x)

        To construct the :math:`\mathbb{L}^2` distance it is used the
        :math:`\mathbb{L}^2` norm wich it is used to compute the distance.

        >>> l2_distance = NormInducedMetric(l2_norm)
        >>> d = l2_distance(fd, fd2)
        >>> float('%.3f'% d)
        0.289

    """

    def __init__(self, norm: Norm[VectorType]):
        self.norm = norm

    def __call__(self, elem1: VectorType, elem2: VectorType) -> NDArrayFloat:
        """Compute the induced norm between two vectors."""
        return self.norm(elem1 - elem2)

    def __repr__(self) -> str:
        return f"{type(self).__name__}(norm={self.norm})"

@multimethod.multidispatch
def pairwise_metric_optimization(
    metric: Any,
    elem1: Any,
    elem2: Optional[Any],
) -> NDArrayFloat:
    r"""
    Optimized computation of a pairwise metric.

    This is a generic function that can be subclassed for different
    combinations of metric and operators in order to provide a more
    efficient implementation for the pairwise metric matrix.
    """
    return NotImplemented

class PairwiseMetric(Generic[MetricElementType]):
    r"""Pairwise metric function.

    Computes a given metric pairwise. The matrix returned by the pairwise
    metric is a matrix with as many rows as observations in the first object
    and as many columns as observations in the second one. Each element
    (i, j) of the matrix is the distance between the ith observation of the
    first object and the jth observation of the second one.

    Args:
        metric: Metric between two elements of a metric
            space.

    """

    def __init__(
        self,
        metric: Metric[MetricElementType],
    ):
        self.metric = metric

    def __call__(
        self,
        elem1: MetricElementType,
        elem2: Optional[MetricElementType] = None,
    ) -> NDArrayFloat:
        """Evaluate the pairwise metric."""
        optimized = pairwise_metric_optimization(self.metric, elem1, elem2)

        return (
            _pairwise_symmetric(self.metric, elem1, elem2)
            if optimized is NotImplemented
            else optimized
        )

    def __repr__(self) -> str:
        return f"{type(self).__name__}(metric={self.metric})"

Original = TypeVar("Original", bound=Vector)
Transformed = TypeVar("Transformed", bound=Vector)

class TransformationMetric(Generic[Original, Transformed], Metric[Original]):
    """
    Compute a distance after transforming the data.

    This is a convenience function to compute a metric after a transformation
    is applied to the data. It can be used, for example, to compute
    Sobolev-like metrics.

    Args:
        e1: First object.
        e2: Second object.

    Returns:
        Distance.

    Examples:
        Compute the L2 distance between the function derivatives.

        >>> import skfda
        >>> from skfda.misc.metrics import l2_distance, TransformationMetric

        >>> x = np.linspace(0, 1, 1001)
        >>> fd = skfda.FDataGrid([x], x)
        >>> fd2 = skfda.FDataGrid([x/2], x)

        >>> dist = TransformationMetric(
        ...     transformation=lambda x: x.derivative(),
        ...     metric=l2_distance,
        ... )
        >>> dist(fd, fd2)
        array([ 0.5])

    """

    def __init__(
        self,
        transformation: Callable[[Original], Transformed],
        metric: Metric[Transformed],
    ):
        self.transformation = transformation
        self.metric = metric

    def __call__(
        self,
        e1: Original,
        e2: Original,
    ) -> NDArrayFloat:
        """Compute the distance."""
        e1_trans = self.transformation(e1)
        e2_trans = self.transformation(e2)

        return self.metric(e1_trans, e2_trans)

    def __repr__(self) -> str:
        return (
            f"{type(self).__name__}()"
        )

@pairwise_metric_optimization.register
def _pairwise_metric_optimization_transformation_distance(
    metric: TransformationMetric[Any, Any],
    e1: T,
    e2: Optional[T],
) -> NDArrayFloat:

    e1_trans = metric.transformation(e1)
    e2_trans = None if e2 is None else metric.transformation(e2)

    pairwise = PairwiseMetric(metric.metric)

    return pairwise(e1_trans, e2_trans)

max_stars_repo_name: jiduque/scikit-fda

max_issues_repo_name: jiduque/scikit-fda

max_forks_repo_name: jiduque/scikit-fda


text:
import math
import sys
import h5py
import numpy
import random
import scipy
from scipy.spatial.distance import cdist

f = h5py.File(sys.argv[1])

# default: relate average distance to distance of 10-th NN
k = 10

if len(sys.argv) > 2:
    k = int(sys.argv[2])

samples = 3000
m = len(f['test'])

distances = numpy.array(f['distances'])

#query_indices = [random.choice(range(len(f['test']))) for _ in range(m)]

queries = numpy.array(f['test'])

dataset = numpy.array(f['train'])

estimates = numpy.zeros(m, dtype=numpy.float)

random_matrix = numpy.array([random.choice(f['train']) for _ in range(samples)])

#for i, v in enumerate(queries):
if 'euclidean' in sys.argv[1]:
    avg = numpy.mean(numpy.transpose(cdist(random_matrix, queries, 'euclidean')), axis=1)
if 'angular' in sys.argv[1]:
    avg = numpy.median(numpy.transpose(cdist(random_matrix, queries, 'cosine')), axis=1)

assert len(avg) == len(queries)

for i in range(len(queries)):
    for j in range(k - 1, 100):
        if distances[i][j] > 1e-6:
            dist = distances[i][j]
            break

    estimates[i] = (avg[i] / dist)

for i,e in enumerate(estimates):
    print(i, e)
# print(sys.argv[1], numpy.median(estimates))

max_stars_repo_name: Cecca/role-of-dimensionality

max_forks_repo_name: Cecca/role-of-dimensionality

reason: import numpy,import scipy,from scipy


text:
import requests
from numpy import NaN
import datetime as datetime
from bs4 import BeautifulSoup
import re
import locale

def getInfo(url):

	item = requests.get(url)
	item = BeautifulSoup(item.text, 'html.parser')

	itemdata = {'itemid'	: NaN,
							'itemname': '',
							'url'			: url,
							'buylimit': NaN,
							'highalch': NaN,
							'lowalch'	: NaN,
							'price'		: NaN,
							'updated'	: datetime.datetime.now,
							'series'	: []
							}

	locale.setlocale( locale.LC_ALL, 'en_US.UTF-8' )

	#grab values from the html
	values = [item.find(id='exchange-itemid').text,
						item.find(class_='gemw-name').text,
						item.find(id='exchange-highalch').text,
						item.find(id='exchange-lowalch').text,
						item.find(class_='gemw-price').text,
						item.find(id='exchange-limit').text
						]

	#try-except blocks to allow default values to override if non-int is passed
	try:
		itemdata['itemid'] = int(values[0])
	except (TypeError, ValueError):
		pass

	if values[1]:
		itemdata['itemname'] = values[1]

	try:
		itemdata['highalch'] = int(locale.atoi(values[2]))
	except (TypeError, ValueError):
		pass

	try:
		itemdata['lowalch'] = int(locale.atoi(values[3]))
	except (TypeError, ValueError):
		pass

	itemdata['price'] = int(locale.atoi(values[4]))

	date = item.find(class_='gemw-updated')['data-date']

	#Fix when able to tell if month is given as Aug or August
	try:
		dt = datetime.datetime.strptime(date, '%d %b %Y, %H:%M (%Z)')
	except ValueError:
		dt = datetime.datetime.strptime(date, '%d %B %Y, %H:%M (%Z)')

	itemdata['updated'] = dt

	try:
		itemdata['buylimit'] = int(locale.atoi(values[5]))
	except (TypeError, ValueError):
		pass

	return itemdata

def getLinks(url, prefix, pattern):

	page = requests.get(url)
	regex = re.compile(pattern)
	pagesoup = BeautifulSoup(page.text, features='lxml')
	links = []
	for link in pagesoup.find_all('a', attrs={'href': regex}):
		links.append(prefix + link.get('href'))
	return links

def getItemPrice(url):

	item = requests.get(url)
	item = BeautifulSoup(item.text, 'html.parser')
	locale.setlocale( locale.LC_ALL, 'en_US.UTF-8' )

	try:
		price = int(locale.atoi(item.find(class_='gemw-price').text))
	except (TypeError, ValueError):
		price = NaN
	return price

max_stars_repo_name: Fradley/rsgetool

max_issues_repo_name: Fradley/rsgetool

max_forks_repo_name: Fradley/rsgetool

reason: from numpy


text:
from __future__ import division
import os
import numpy as np
import h5py
import keras.backend as K
from keras.callbacks import Callback

class ReplayData(Callback):
    """Creates an instance of a Keras Callback to collect information
    about the training process, mostly the weights, and store them
    in a group of a HDF5 file, together with inputs and targets passed
    as arguments.

    It also saves the Keras model itself, as an independent HDF5 file,
    named after the group_name with '_model' appended to it.

    Parameters
    ----------
    inputs: ndarray
        An array with the inputs to be used during the training of the
        model. These inputs are going to be later used to generate
        statistics and plots.
    targets: ndarray
        An array with the targets to be used during the training of the
        model. These targets are going to be later used to generate
        statistics and plots.
    filename: String
        HDF5 filename to be used to store the collected information.
        It can be an existing file, as the data will be appended to it.
    group_name: String
        Group inside the HDF5 file where the information is to be
        saved. If the informed group name already exists, it will throw
        an exception.
    model: Keras Model, optional
        If provided, it will set the model directly to the callback
        instance and execute `on_train_begin` method to initialize
        all variables and create the corresponding group in the HDF5
        file.
        This is intended to be used for analyzing the initial conditions
        of the model without ever calling its `fit` function, where
        the callback is usually called.
    """
    def __init__(self, inputs, targets, filename, group_name, model=None):
        super(ReplayData, self).__init__()
        self.handler = h5py.File('{}'.format(filename), 'a')
        self.inputs = inputs
        self.targets = targets.reshape(len(targets), -1)
        self.filepath = os.path.split(filename)[0]
        self.filename = filename
        self.group = None
        self.group_name = group_name
        self.current_epoch = -1
        self.n_epochs = 0
        if model is not None:
            self.set_model(model)
            self.set_params({
                'epochs': 0,
                'samples': len(self.inputs),
                'batch_size': len(self.inputs),
            })
            self.group_name = group_name + '_init'
            self.on_train_begin()
            self.group_name = group_name
        return

    def _append_weights(self):
        """Appends the weights to the corresponding layer and weight sequence inside the group, for the current epoch.
        """
        for i, layer in enumerate(self.model.layers):
            layer_weights = layer.get_weights()
            for j, weights in enumerate(layer_weights):
                self.group['layer{}'.format(i)]['weights{}'.format(j)][self.current_epoch + 1] = weights

    def get_lr(self):
        optimizer = self.model.optimizer
        return K.function(inputs=[],
                          outputs=[optimizer.lr *
                                   (1. / (1. + optimizer.decay * K.cast(optimizer.iterations,
                                                                        K.dtype(optimizer.decay))))])(inputs=[])[0]

    def on_train_begin(self, logs={}):
        self.model.save(os.path.join(self.filepath, '{}_model.h5'.format(self.group_name)))
        self.n_epochs = self.params['epochs']

        self.group = self.handler.create_group(self.group_name)
        self.group.attrs['samples'] = self.params['samples']
        self.group.attrs['batch_size'] = self.params['batch_size']
        self.group.attrs['n_batches'] = (self.params['samples'] + self.params['batch_size'] - 1) // \
                                        self.params['batch_size']
        self.group.attrs['n_epochs'] = self.n_epochs
        self.group.attrs['n_layers'] = len(self.model.layers)
        try:
            # Python 2
            self.group.attrs['activation_functions'] = [layer.activation.func_name
                                                        if hasattr(layer, 'activation')
                                                        else ''
                                                        for layer in self.model.layers]
        except AttributeError:
            # Python 3
            self.group.attrs['activation_functions'] = [np.string_(layer.activation.__name__)
                                                        if hasattr(layer, 'activation')
                                                        else np.string_('')
                                                        for layer in self.model.layers]
        self.group.create_dataset('inputs', data=self.inputs)
        self.group.create_dataset('targets', data=self.targets)

        self.group.create_dataset('loss', shape=(self.n_epochs,), dtype='f')
        for metric in self.model.metrics:
            metric_name = metric
            if callable(metric):
                metric_name = metric.name if hasattr(metric, 'name') else metric.__name__
            self.group.create_dataset(metric_name, shape=(self.n_epochs,), dtype='f')

        self.group.create_dataset('lr', shape=(self.n_epochs,), dtype='f')

        for i, layer in enumerate(self.model.layers):
            layer_grp = self.group.create_group('layer{}'.format(i))
            layer_weights = layer.get_weights()
            for j, weights in enumerate(layer_weights):
                layer_grp.create_dataset('weights{}'.format(j),
                                         shape=(self.n_epochs + 1, ) + weights.shape,
                                         dtype='f')
        self._append_weights()
        return

    def on_train_end(self, logs={}):
        self.handler.close()
        return

    def on_epoch_begin(self, epoch, logs={}):
        self.current_epoch = epoch
        self.group['lr'][epoch] = self.get_lr()
        return

    def on_epoch_end(self, epoch, logs={}):
        self._append_weights()
        self.group['loss'][epoch] = logs.get('loss')
        for metric in self.model.metrics:
            metric_name = metric
            if callable(metric):
                metric_name = metric.name if hasattr(metric, 'name') else metric.__name__
            self.group[metric_name][epoch] = logs.get(metric_name, np.nan)
        return

max_stars_repo_name: dackerer/deepreplay

max_issues_repo_name: dvirginz/deepreplay

max_forks_repo_name: dvirginz/deepreplay


text:
########################################################################
# Copyright 2021, UChicago Argonne, LLC
#
# Licensed under the BSD-3 License (the "License"); you may not use
# this file except in compliance with the License. You may obtain a
# copy of the License at
#
#     https://opensource.org/licenses/BSD-3-Clause
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or
# implied. See the License for the specific language governing
# permissions and limitations under the License.
########################################################################
"""
date: 2020-06-09
author: matz
Upgraded Cheng-Todreas Detailed correlations (2018)
"""
########################################################################
import numpy as np
from . import friction_ctd as ctd

# Application ranges of friction factor correlations
applicability = {}
applicability['P/D'] = np.array([1.0, 1.42])
applicability['H/D'] = np.array([8.0, 52.0])
applicability['Nr'] = np.array([7, 217])
applicability['regime'] = ['turbulent', 'transition', 'laminar']
applicability['Re'] = np.array([50, 1e6])
applicability['bare rod'] = True

########################################################################
# CONSTANTS - based only on assembly geometry
########################################################################

def calc_constants(asm_obj):
    """Calculate and store constants for friction factor calculation

    Parameters
    ----------
    asm_obj : DASSH Assembly object
        Contains geometric parameters on which constants are based

    Returns
    -------
    dict
        Constants used in friction factor calculation:
        1. Reynolds number flow regime bounds {list}
            (laminar-transition-turbulent)
        2. Subchannel friction factor constants {dict}
            (turbulent, laminar; one for each type of subchannel)
        3. Bundle friction factor constants {dict}
            (turbulent, laminar)

    """
    c = {}
    c['Re_bnds'] = calculate_Re_bounds(asm_obj)
    c['Cf_sc'] = calculate_subchannel_friction_factor_const(asm_obj)
    c['Cf_b'] = calculate_bundle_friction_factor_const(asm_obj)
    return c

########################################################################
# REYNOLDS NUMBER BOUNDARIES: laminar-transition, transition-turbulent
########################################################################

def calculate_Re_bounds(asm_obj):
    """Calculate the Reynolds numbers for the boundaries of the
    laminar, transition, and turbulent flow regimes

    Parameters
    ----------
    asm_obj : DASSH assembly object
        Required attributes: pin diameter (m) and pin pitch (m)

    Returns
    -------
    tuple
        Reynolds numbers (float) for the boundaries between the
        laminar-transition and transition-turbulent flow regimes

    """
    re_bt = 1e4 * 10**(0.7 * (asm_obj.pin_pitch
                              / asm_obj.pin_diameter - 1.0))
    re_bl = 320.0 * 10**(asm_obj.pin_pitch
                         / asm_obj.pin_diameter - 1.0)
    return (re_bl, re_bt)

########################################################################
# SUBCHANNEL FRICTION FACTOR CONSTANTS
########################################################################

def calculate_subchannel_friction_factor_const(asm):
    """Calculate the Upgraded Cheng-Todreas Detailed correlation
    subchannel friction factor constants (Cf)

    Parameters
    ----------
    asm : DASSH Assembly object
        Contains the assembly geometry data

    Returns
    -------
    dict
        Friction factor constants (float) for different coolant
        subchannels in different flow regimes (keys: ["laminar",
        "turbulent"])
    """
    # Get the Cheng-Todreas friction constant polynomial factors
    # and other constants required for the correlation
    wire_sweep = calculate_wire_sweeping_const(asm.wire_pitch
                                               / asm.pin_diameter)
    wire_drag = calculate_wire_drag_const(asm.wire_diameter
                                          / asm.pin_diameter,
                                          asm.wire_pitch
                                          / asm.pin_diameter)
    # Calculate friction factor constants
    return ctd._calc_sc_ff_const(asm, wire_drag, wire_sweep)

def calculate_wire_sweeping_const(h_over_d):
    """Calculate the wire-sweeping constant for the Upgraded
    Cheng-Todreas friction factor constant calculation"""
    ws = {}
    if h_over_d == 0.0:
        ws['turbulent'] = 0.0
        ws['laminar'] = 0.0
    else:
        ws['turbulent'] = -11.0 * np.log10(h_over_d) + 19.0
        ws['laminar'] = ws['turbulent']
    return ws

def calculate_wire_drag_const(wd_over_d, h_over_d):
    """Calculate the wire-drag constant for the Upgraded
    Cheng-Todreas friction factor constant calculation"""
    wd = {}
    if wd_over_d == 0.0:
        wd['turbulent'] = 0.0
        wd['laminar'] = 0.0
    else:
        wd['turbulent'] = ((19.56 - 98.71 * wd_over_d + 303.47 * wd_over_d**2)
                           / h_over_d**0.541)
        wd['laminar'] = 1.4 * wd['turbulent']
    return wd

########################################################################
# BUNDLE FRICTION FACTOR CONSTANT
########################################################################

def calculate_bundle_friction_factor_const(asm_obj):
    """Calculate the bundle (assembly avg) friction factor constant

    Parameters
    ----------
    asm_obj : DASSH Assembly object
        Contains geometric and flow parameters
    subchannel_cf : dict
        dict (keys: ['laminar', 'turbulent']) of numpy.ndarray
        (length = 3), each containing friction factor constants
        for each coolant subchannel

    Notes
    -----
    This is a wrapper function used to apply UCTD (2018) parameters
    to the bundle average friction factor constant calculation

    """
    try:
        cf_sc = asm_obj.corr_constants['ff']['Cf_sc']
    except (KeyError, AttributeError):
        cf_sc = calculate_subchannel_friction_factor_const(asm_obj)

    return ctd._calc_cfb(asm_obj, cf_sc)

########################################################################
# BUNDLE FRICTION FACTOR
########################################################################

def calculate_bundle_friction_factor(asm_obj):
    """Calculate the bundle-average friction factor

    Parameters
    ----------
    asm_obj : DASSH Assembly object
        Contains the assembly geometric details and bundle Re

    Returns
    -------
    float
        Bundle-average friction factor based on assembly geometry
        and flow regime

    """
    try:
        cfb = asm_obj.corr_constants['ff']['Cf_b']
    except (KeyError, AttributeError):
        cfb = calculate_bundle_friction_factor_const(asm_obj)

    f = {}
    # Calculate friction factor for laminar and turbulent regimes
    for r in cfb.keys():
        f[r] = cfb[r] / asm_obj.coolant_int_params['Re']**ctd._m[r]
    # If transition region, combine laminar and turbulent friction
    # factors using intermittency fxn; otherwise, return value
    try:
        Re_bl, Re_bt = asm_obj.corr_constants['ff']['Re_bnds']
    except (KeyError, AttributeError):
        Re_bl, Re_bt = calculate_Re_bounds(asm_obj)

    if asm_obj.coolant_int_params['Re'] <= Re_bl:
        return f['laminar']
    elif asm_obj.coolant_int_params['Re'] >= Re_bt:
        return f['turbulent']
    else:
        # Transition regime intermittency factor
        x = ctd.calc_intermittency_factor(asm_obj, Re_bt, Re_bl)
        # Different correlation for transition region than O.G. CTD
        return (f['laminar'] * (1 - x)**(1 / 3.0) * (1 - x**7.0)
                + f['turbulent'] * x**(1 / 3.0))

max_stars_repo_name: khurrumsaleem/dassh

max_issues_repo_name: khurrumsaleem/dassh

max_forks_repo_name: khurrumsaleem/dassh