text:
import numpy as np
from modules.containers import VoigtLine
from modules.redshifter import RedshiftMeasurer
from modules.spectrumgen_voigts import VoigtLinesGenerator
from modules.rvsimulator_singletimeorbit import SingleTimeOrbitRVSimulator
from modules.s1reader import S1Reader
from modules.settings import Settings
from modules.helpers import ensurePathExists, savedatafile, truncate
from modules.s1spectrafetcher import S1SpectraFetcher

class RVSpectraSimulator():
	def __init__(self, Mstar=1.0, Mplanet=1.0, a=1.0, e=0.0, w=0.0, v0=0.0, SNR=0.0, R=1e5, times=np.arange(0.0, 365.25, 0.1)):
		self.G = 8.88677e-10 # AU^3 Mearth^-1 day^-2
		self.SunInEarthMasses = 332978.9 # 1 solar mass in Earth masses

		self.Mstar = Mstar # solar masses
		self.Mplanet = Mplanet # Earth masses
		self.a = a # AU
		self.e = e # 0..1
		self.w = w # argument of periastron, radians
		self.v0 = v0 # barycenter velocity (constant component)

		self.SNR = SNR # set to 0 to eliminate all noise
		self.R = R # spectral resolution Lambda/dLambda
		self.times = times

		self.simulator = SingleTimeOrbitRVSimulator(times, Mstar, Mplanet, a, e, w, v0)
		self.redshiftMeasurer = RedshiftMeasurer()

		self.lines = np.asarray([])
		self.continuum = 1.0
		self.absorption = True
		self.margin = 5 # angstroems

		self.maxshift = 1e3 # m/s
		self.sfactor = 300 # for spline interpolation of the line

		self.settings = Settings()

	def _getLinesProperties(self):
		if len(self.lines) == 0:
			raise Exception('Provide at least one item of type <VoigtLine> for the spectrum in <RVSpectraSimulator.lines[]>')

		centers = [VoigtLine.center for VoigtLine in self.lines]
		fluxes = [VoigtLine.flux for VoigtLine in self.lines]
		g_fwhms = [VoigtLine.g_fwhm for VoigtLine in self.lines]
		l_fwhms = [VoigtLine.l_fwhm for VoigtLine in self.lines]

		return centers, fluxes, g_fwhms, l_fwhms

	def _getModelSpectrum(self):
		centers, fluxes, g_fwhms, l_fwhms = self._getLinesProperties()
		model = VoigtLinesGenerator(centers=centers, margin=self.margin, rv=0, g_fwhm=g_fwhms, l_fwhm=l_fwhms, R=self.R, fluxes=fluxes, cont=self.continuum)

		return model.getSpectrum(absorption=self.absorption)

	def addLine(self, voigt):
		self.lines = np.append(self.lines, voigt)

	def removeAllLines(self):
		self.lines = np.asarray([])

	def period(self):
		return np.power(4.0*np.power(np.pi, 2)*np.power(self.a, 3) / (self.G*(self.Mstar * self.SunInEarthMasses + self.Mplanet)), 0.5)

	def getRVCurve(self, times=[], dt=None):
		if dt != None:
			times = np.arange(min(self.times), max(self.times), dt)

		if len(times) > 0:
			simulator = SingleTimeOrbitRVSimulator(times, self.Mstar, self.Mplanet, self.a, self.e, self.w, self.v0)
			return simulator.getRVCurve()

		return self.simulator.getRVCurve()

	def getRVData(self, progressbar=False, template='average', outfunction=None, times_full=None, cycle=True):
		ts0, rvs0 = self.getRVCurve()
		fetcher = S1SpectraFetcher()
		time_from = self.settings.time_from

		wl_model, flux_model = [], []
		lines_wls = []
		lines_fluxes = []

		if progressbar and outfunction != None:
			outfunction( 'producing profiles ...' )
		# get lines first

		points = []
		if progressbar:
			from progressbar import ProgressBar,SimpleProgress,ETA
			pbar = ProgressBar(widgets=['                        ', SimpleProgress(), ', ', ETA()])
			points = pbar(range(0,len(rvs0)))
		else:
			points = range(0,len(rvs0))

		for rvi in points:
			# stellar rotation phase
			start_rphase = self.settings.start_rphase
			rotation_period = self.settings.rotation_period
			rphase = start_rphase + (ts0[rvi] - time_from) / rotation_period

			# magnetic cycle phase
			start_mphase = self.settings.start_mphase
			magnetic_period = self.settings.magnetic_period
			mphase = start_mphase + (ts0[rvi] - time_from) / magnetic_period

			# get profile
			wls, fluxes = fetcher.getProfile_old(rphase, mphase, rv=rvs0[rvi], timepoint=ts0[rvi])
			#wls, fluxes = fetcher.getProfile(time=ts0[rvi], rphase=rphase, rv=rvs0[rvi])

			if self.SNR > 0.0:
				# add Gaussian noise with sigma = 1.0/SNR
				noise = np.random.normal(0, 1.0/self.SNR, (len(fluxes)))
				fluxes = np.asarray([x+y for x, y in zip(fluxes, noise)])

			# save profile to file
			if len(self.settings.profiles_path) > 0:
				path = self.settings.profiles_path.replace(r'{N}', str(rvi)).replace(r'{RV}', str(truncate(rvs0[rvi]))).replace(r'{time}', str(truncate(ts0[rvi])))
				ensurePathExists(path)
				lsd_rvs = wls#reader.wlToRv(phase=rphase, wls=wls)
				savedatafile(path, \
					(wls, lsd_rvs, fluxes), ('wl[A]', 'rv_shift[m/s]', 'flux'), \
					comment='RV = ' + str(rvs0[rvi]) + ' [m/s],\tt = ' + str(ts0[rvi]) + ' [days],\tstellar rotation phase: ' + str(rphase))

			lines_wls.append(wls)
			lines_fluxes.append(fluxes)

		self.R = fetcher.R

		if progressbar and outfunction != None:
			outfunction( 'computing template ...' )

		if template == 'average':
			# average lines and use it as a model
			for rvi in range(0,len(rvs0)):
				if rvi == 0:
					for wli in range(0,len(lines_wls[rvi])):
						wl_model.append(lines_wls[rvi][wli])
						flux_model.append(lines_fluxes[rvi][wli])
				else:
					for wli in range(0,len(lines_wls[rvi])):
						flux_model[wli] += lines_fluxes[rvi][wli]

				lines_wls[rvi] = np.asarray(lines_wls[rvi])
				lines_fluxes[rvi] = np.asarray(lines_fluxes[rvi])
			n = len(rvs0)
			for i in range(0, len(flux_model)):
				flux_model[i] /= n

			wl_model = np.asarray(wl_model)
			flux_model = np.asarray(flux_model)

		elif template == 'model':
			wl_model, flux_model = self._getModelSpectrum()
		else:
			raise Exception('Unknown template selected in RVSpectraSimulator.getRVData(). Available options: "model" and "average".')

		if progressbar and outfunction != None:
			outfunction( 'measuring RVs ...' )

		points = []
		if progressbar:
			from progressbar import ProgressBar,SimpleProgress,ETA
			pbar = ProgressBar(widgets=['                        ', SimpleProgress(), ', ', ETA()])
			points = pbar(rvs0)
		else:
			points = rvs0

		i = 0
		rvdata = []
		for rv in points:
			# measure redshift
			val = self.redshiftMeasurer.getRV(lines_wls[i], lines_fluxes[i], wl_model, flux_model, maxshift=self.maxshift, R=self.R, fit_type='linespline', sfactor=self.sfactor)
			rvdata.append(val)
			i += 1

		return ts0, np.asarray(rvdata), ts0, rvs0

max_stars_repo_name: timberhill/radiant

max_issues_repo_name: timberhill/radiant

max_forks_repo_name: timberhill/radiant


text:
# Testing raytracing functions against lenstools
import tensorflow as tf
import numpy as np
from numpy.testing import assert_allclose
from scipy.interpolate import InterpolatedUnivariateSpline as iuspline
from scipy.ndimage import fourier_gaussian
import flowpm
import flowpm.constants as constants
import lenstools as lt
import bigfile
import os
import astropy.units as u

data_path = os.path.join(os.path.dirname(os.path.abspath(__file__)), '..')

np.random.seed(0)

bs = 200
nc = 64
plane_resolution = 64
field = 5
npix = 64

class FlowPMSnapshot(lt.simulations.nbody.NbodySnapshot):
  """
    A class that handles FlowPM simulation snapshots for lenstools
    """

  _header_keys = [
      'masses', 'num_particles_file', 'num_particles_total', 'box_size',
      'num_files', 'Om0', 'Ode0', 'h'
  ]

  ############################
  #Open the file with bigfile#
  ############################

  @classmethod
  def open(cls, filename, pool=None, header_kwargs=dict(), **kwargs):

    if bigfile is None:
      raise ImportError("bigfile must be installed!")

    fp = bigfile.BigFile(cls.buildFilename(filename, pool, **kwargs))
    return cls(fp, pool, header_kwargs=header_kwargs)

  ###################################################################################
  ######################Abstract method implementation###############################
  ###################################################################################

  @classmethod
  def buildFilename(cls, root, pool):
    return root

  @classmethod
  def int2root(cls, name, n):
    return name

  def getHeader(self):

    #Initialize header
    header = dict()
    bf_header = self.fp["Header"].attrs

    ###############################################
    #Translate FlowPM header into lenstools header#
    ###############################################

    #Number of particles/files
    header["num_particles_file"] = bf_header["NC"][0]**3
    header["num_particles_total"] = header["num_particles_file"]
    header["num_files"] = 1

    #Cosmology
    header["Om0"] = bf_header["OmegaCDM"][0] + bf_header["OmegaB"][0]
    header["Ode0"] = 1. - header["Om0"]
    header["w0"] = -1.
    header["wa"] = 0.
    header["h"] = bf_header["h"][0]
    header["redshift"] = 1 / (bf_header["Time"][0]) - 1
    header["comoving_distance"] = bf_header["comoving_distance"][0] * 1.0e3
    header["scale_factor"] = bf_header["Time"][0]

    #Box size in kpc/h
    header["box_size"] = bf_header["BoxSize"][0] * 1.0e3
    header["box_size_mpch"] = bf_header["BoxSize"][0]
    #Plane Resolution
    header["nc"] = bf_header["NC"][0]
    #Masses
    header["masses"] = np.array(
        [0., bf_header["M0"][0] * header["h"], 0., 0., 0., 0.])
    #################

    return header

  def setLimits(self):

    if self.pool is None:
      self._first = None
      self._last = None
    else:

      #Divide equally between tasks
      Nt, Np = self.pool.size + 1, bigfile.BigData(self.fp).size
      part_per_task = Np // Nt
      self._first = part_per_task * self.pool.rank
      self._last = part_per_task * (self.pool.rank + 1)

      #Add the remainder to the last task
      if (Np % Nt) and (self.pool.rank == Nt - 1):
        self._last += Np % Nt

  def getPositions(self, first=None, last=None, save=True):

    #Get data pointer
    data = self.fp

    #Read in positions in Mpc/h
    if (first is None) or (last is None):
      positions = (data["0/Position"][:] +
                   np.array([
                       0.5 / self.header["nc"] * self.header["box_size_mpch"],
                       0.5 / self.header["nc"] * self.header["box_size_mpch"], 0
                   ],
                            dtype=np.float32)) * self.Mpc_over_h
    else:
      positions = data["0/Position"][first:last] * self.Mpc_over_h

    #Enforce periodic boundary conditions
    for n in (0, 1):
      positions[:, n][positions[:, n] < 0] += self.header["box_size"]
      positions[:, n][
          positions[:, n] > self.header["box_size"]] -= self.header["box_size"]

    #Maybe save
    if save:
      self.positions = positions

    #Initialize useless attributes to None
    self.weights = None
    self.virial_radius = None
    self.concentration = None

    #Return
    return positions

def test_density_plane(return_results=False):
  """ Tests cutting density planes from snapshots against lenstools
  """
  klin = np.loadtxt(data_path + '/flowpm/data/Planck15_a1p00.txt').T[0]
  plin = np.loadtxt(data_path + '/flowpm/data/Planck15_a1p00.txt').T[1]
  ipklin = iuspline(klin, plin)

  cosmo = flowpm.cosmology.Planck15()

  a0 = 0.9
  r0 = flowpm.background.rad_comoving_distance(cosmo, a0)

  # Create a state vector
  initial_conditions = flowpm.linear_field(nc, bs, ipklin, batch_size=2)
  state = flowpm.lpt_init(cosmo, initial_conditions, a0)

  # Export the snapshot
  flowpm.io.save_state(
      cosmo,
      state,
      a0, [nc, nc, nc], [bs, bs, bs],
      'snapshot_density_testing',
      attrs={'comoving_distance': r0})

  # Reload the snapshot with lenstools
  snapshot = FlowPMSnapshot.open('snapshot_density_testing')

  # Cut a lensplane in the middle of the volume
  lt_plane, resolution, NumPart = snapshot.cutPlaneGaussianGrid(
      normal=2,
      plane_resolution=plane_resolution,
      center=(bs / 2) * snapshot.Mpc_over_h,
      thickness=(bs / 4) * snapshot.Mpc_over_h,
      left_corner=np.zeros(3) * snapshot.Mpc_over_h,
      smooth=None,
      kind='density')

  # Cut the same lensplane with flowpm
  fpm_plane = flowpm.raytracing.density_plane(
      state, nc, center=nc / 2, width=nc / 4, plane_resolution=plane_resolution)

  # Apply additional normalization terms to match lenstools definitions
  constant_factor = 3 / 2 * cosmo.Omega_m * (constants.H0 / constants.c)**2
  density_normalization = bs / 4 * r0 / a0
  fpm_plane = fpm_plane * density_normalization * constant_factor

  # Checking first the mean value, which accounts for any normalization
  # issues
  assert_allclose(np.mean(fpm_plane[0]), np.mean(lt_plane), rtol=1e-5)

  # To check pixelwise difference, we need to do some smoothing as lenstools and
  # flowpm use different painting kernels
  smooth_lt_plane = np.fft.ifft2(fourier_gaussian(np.fft.fft2(lt_plane),
                                                  3)).real
  smooth_fpm_plane = np.fft.ifft2(
      fourier_gaussian(np.fft.fft2(fpm_plane[0]), 3)).real

  assert_allclose(smooth_fpm_plane, smooth_lt_plane, rtol=2e-2)

  if return_results:
    return fpm_plane, lt_plane, smooth_fpm_plane, smooth_lt_plane

def test_convergence_Born(return_results=False):
  """ This function tests that given a set of density planes,
  both lenstools and flowpm recover the same convergence maps in
  angular coordinates.
  """
  klin = np.loadtxt(data_path + '/flowpm/data/Planck15_a1p00.txt').T[0]
  plin = np.loadtxt(data_path + '/flowpm/data/Planck15_a1p00.txt').T[1]
  ipklin = iuspline(klin, plin)

  cosmo = flowpm.cosmology.Planck15()

  a0 = 0.9

  # Create a state vector
  initial_conditions = flowpm.linear_field([nc, nc, 10 * nc], [bs, bs, 10 * bs],
                                           ipklin,
                                           batch_size=2)
  state = flowpm.lpt_init(cosmo, initial_conditions, a0)

  r = tf.linspace(0., 10 * bs, 11)
  r_center = 0.5 * (r[1:] + r[:-1])
  a_center = flowpm.background.a_of_chi(cosmo, r_center)

  constant_factor = 3 / 2 * cosmo.Omega_m * (constants.H0 / constants.c)**2

  # To make it convenient to access simulation properties in lenstools
  # let's quicly export and reload the sim
  # TODO: remove the need for this!
  flowpm.io.save_state(
      cosmo,
      state,
      a0, [nc, nc, 10 * nc], [bs, bs, 10 * bs],
      'snapshot_born_testing',
      attrs={'comoving_distance': r_center[0]})

  # Reload the snapshot with lenstools
  snapshot = FlowPMSnapshot.open('snapshot_born_testing')

  # Get some density planes and create lenstool tracer
  lensplanes = []
  tracer = lt.simulations.RayTracer(lens_type=lt.simulations.DensityPlane)
  for i in range(len(r_center)):
    plane = flowpm.raytracing.density_plane(
        state, [nc, nc, 10 * nc],
        r_center[i] / bs * nc,
        width=nc,
        plane_resolution=plane_resolution)
    r, a, p = r_center[i], a_center[i], plane[0]
    lensplanes.append((r, a, plane))

    density_normalization = bs * r / a

    # We upsample the lensplanes before giving them to lenstools because
    # lentools is using a weird kind of interpolation when converting from
    # comoving coordinates to angular coords. with a larger
    p = tf.image.resize(
        tf.reshape(p, [1, plane_resolution, plane_resolution, 1]), [2048, 2048])

    p = (p[0, :, :, 0] * constant_factor * density_normalization).numpy()
    p = p - np.mean(p)
    lt_plane = lt.simulations.DensityPlane(
        p,
        angle=snapshot.header["box_size"],
        redshift=1 / a - 1,
        cosmology=snapshot.cosmology)
    tracer.addLens(lt_plane)

  # Adding dummy lensplane at the end
  tracer.addLens(
      lt.simulations.DensityPlane(
          np.zeros((2048, 2048)),
          angle=snapshot.header["box_size"],
          redshift=0.99,
          cosmology=snapshot.cosmology))
  tracer.addLens(
      lt.simulations.DensityPlane(
          np.zeros((2048, 2048)),
          angle=snapshot.header["box_size"],
          redshift=2,
          cosmology=snapshot.cosmology))
  tracer.reorderLenses()

  # Create an array of coordinates at which to retrieve the convernge maps
  xgrid, ygrid = np.meshgrid(
      np.linspace(0, field, npix, endpoint=False),  # range of X coordinates
      np.linspace(0, field, npix, endpoint=False))  # range of Y coordinates

  coords = np.stack([xgrid, ygrid], axis=0) * u.deg
  c = coords.reshape([2, -1]).T

  # Compute convergence map with lenstool
  lt_map = tracer.convergenceBorn(coords, z=1.0)

  # Compute convergemce map with flowpm
  fpm_map = flowpm.raytracing.convergenceBorn(
      cosmo, lensplanes, bs / nc, bs, c.to(u.rad), z_source=tf.ones([1]))

  # Comparing the final maps
  assert_allclose(
      lt_map, fpm_map[0].numpy().reshape([npix, npix, -1])[:, :, -1], atol=5e-4)

  if return_results:
    return lt_map, fpm_map

ext: py

max_stars_repo_name: LBJ-Wade/flowpm

max_issues_repo_name: LBJ-Wade/flowpm

max_forks_repo_name: LBJ-Wade/flowpm

reason: import numpy,from numpy,from scipy,import astropy


text:
"""
extra the last layer embedding of inception3

refer: https://github.com/pytorch/examples/blob/master/imagenet/main.py
"""

import argparse
from tqdm import tqdm
import numpy as np
import json
import datetime
import sys
import os
from scipy.io import loadmat
import pickle

import torch.nn.functional as F
import torch.nn as nn
import torch
import yaml
import torchvision.transforms as transforms
import networkx as nx
#import torchvision.datasets as datasets
from torch.autograd import Variable

from tensorboardX import SummaryWriter

import utils
from models.S2ILinearModel import S2ILayer2Model_IAS
TMP=10

class train(object):
    def __init__(self,device,cfg_file,IAS=False):
        self.IAS = IAS
        self.device = device
        with open(cfg_file,'r') as f:
           self.cfg = yaml.safe_load(f)
        utils.init_seeds(self.cfg['model_hyp']['random_seed'])

        attSplit = loadmat('./dataset/xlsa/'+self.cfg['dataset']['name']+'/att_splits.mat')
        res101 = loadmat('./dataset/xlsa/'+self.cfg['dataset']['name']+'/res101.mat')
        labels = res101['labels'].astype(int).squeeze() - 1
        seen_dataLoc = attSplit['test_seen_loc'].squeeze() - 1
        unseen_dataLoc = attSplit['test_unseen_loc'].squeeze() -1

        seen_labels = labels[seen_dataLoc]
        unseen_labels = labels[unseen_dataLoc]
        self.seen_labels = np.unique(seen_labels)
        self.unseen_labels = np.unique(unseen_labels)

        self.clsname = [ attSplit['allclasses_names'][i][0][0] for i in range(len(attSplit['allclasses_names']))]
        att_matrix = np.transpose(attSplit['att'])
        self.cfg['model_hyp']['att_feats'] = att_matrix.shape[1]

        self.attMatrix = att_matrix.copy()
        self.attMatrix[:len(self.seen_labels)] = att_matrix[self.seen_labels]
        self.attMatrix[len(self.seen_labels):] = att_matrix[self.unseen_labels]

        self.attMatrix = torch.FloatTensor(self.attMatrix).to(self.device)

    def createLinearModel(self,):

        self.model1 = S2ILayer2Model_IAS(att_dims=self.cfg['model_hyp']['att_feats'],img_dims=self.cfg['model_hyp']['img_feats'],tmp=5).to(self.device)

    def apply_classification_weights(self,features, cls_weights,norm=False):

        features = F.normalize(features,dim=-1)
        cls_weights = F.normalize(cls_weights, p=2, dim=-1, eps=1e-12)

        cls_scores = self.scale_cls * (torch.matmul(cls_weights,features.t()))
        cls_scores = cls_scores.permute(0,2,1)
        cls_scores = torch.diagonal(cls_scores,offset=0,dim1=0,dim2=1)
        cls_scores = cls_scores.t()
        if norm:
            return cls_scores,cls_scores/scale_cls
        else:
                return cls_scores

    def IASatt(self,features,AttM,tmp):
        cls_num = len(AttM)
        attdims = AttM.shape[1]
        atten = torch.mm(features,self.w_IAS)+self.b_IAS
        atten = F.relu(atten)
        atten = F.softmax(atten/tmp,dim=1).reshape(-1,1,attdims)
        atten = atten + torch.ones_like(atten)
        atten = atten.repeat(1,cls_num,1)
        AttM = AttM.unsqueeze(0)
        AttM = AttM.repeat(len(features),1,1)
        AttM = atten*AttM
        return AttM

    def forward(self,att,features,tmp):

        features = features.squeeze()
        att = self.IASatt(features,att,tmp)

        a1 = F.relu(torch.matmul(att, self.w1) + self.b1)
        a2 = F.relu(torch.matmul(a1, self.w2) + self.b2)
        return a2

    def loadDataFromCustom(self,dataType):
        """
            Load Data from Xlsa:
            ------

        """
        print('load from customed')

        with open('./dataset/finetune/'+self.cfg['dataset']['name']+'/dvbeExtracted.pkl','rb') as f:
            ExtractFeature = pickle.load(f)
        imgfs = np.array(ExtractFeature[dataType]['features'])
        labels = np.array(ExtractFeature[dataType]['labels'])

        print(imgfs.shape)
        print(labels.shape)
        select_id, idx = np.unique(labels,return_inverse=True)
        if dataType.lower() == "test_unseen":
            labels = idx + len(self.seen_labels)

        else:
            labels = idx

        imgfs = torch.FloatTensor(imgfs.reshape(-1,2048))
        labels = torch.LongTensor(labels)
        dataset = torch.utils.data.TensorDataset(imgfs,labels)
        dataloader = torch.utils.data.DataLoader(dataset,batch_size=self.cfg['dataset_hyp']['batch_size'],shuffle=False)
        return dataloader

    def loadChk(self,chkFile1,chkFile2):
        """
            if resume training
            args:
                chkFile
            return:
                bestAcc in chkFile
        """

        print("=> load checkpoint '{}'".format(chkFile1))
        chk = torch.load(chkFile1)
        self.cfg['model_hyp']['start_epoch'] = chk['epoch']
        print(chk['model'].keys())
        self.model1.load_state_dict(chk['model'])

        print("=> load checkpoint '{}'".format(chkFile2))

        model = torch.load(chkFile2)
        self.w1 = Variable(model['w1'], requires_grad=False)
        self.b1 = Variable(model['b1'], requires_grad=False)
        self.w2 = Variable(model['w2'], requires_grad=False)
        self.b2 = Variable(model['b2'], requires_grad=False)
        self.w_IAS = Variable(model['w_IAS'],requires_grad=False)
        self.b_IAS = Variable(model['b_IAS'],requires_grad=False)
        self.scale_cls = Variable(model['scale_cls'], requires_grad=False)
        self.bias = Variable(model['bias'], requires_grad=False)
        return

    def evaluationSplit(self,dataloader,dataType='seen',accType='allCls',imgfType='customed',recordFile=None,seenOrunseenCls=None):
        """
        evaluation model on input dataloader
        args:
            dataloader
            dataType: 'seen'   : count acc, CEloss, MSEloss
                      'unseen' : count acc, CEloss(1), MSEloss
            accType:  'allCls' : return acc mean of all cls
                      'allImg' : return acc mean of all images
            recordFile: None   : do not record
                        'print'： print eval
                       filepath: record in this file
        return:
            accs: (int) mean of all cls or all images
                        decide by accType
            accs_cls: (dict) mean acc according to cls
            generalAccs_cls: (dict) mean generalAccs 2 cls,
            splitacc: datasplit bw seen and unseen
            cos_cls:(dict){cls:[sample_cos]} sample cosimilarity
        """

        total_spp_loss = []
        generalAccs = []
        accs = []
        generalAccs_cls = {}
        accs_cls = {}
        cos_cls = {}
        spp_loss_cls = {}
        splitaccs = []
        TF_l = []
        cos = []

        with torch.no_grad():
            for images,targets in tqdm(dataloader):
                images = images.to(self.device)
                targets = targets.to(self.device)

                if self.cfg['dataset_hyp']['imgfType'] == 'customed':
                    images = self.caffeRes101(images).squeeze()

                batch_visual_norm = F.normalize(images, p=2, dim=images.dim()-1, eps=1e-12)
                generalpreds1,generalpreds1_n = self.model1(self.attMatrix,images,TrainOrTest="Test")
                batch_weights = self.forward(self.attMatrix[len(self.seen_labels):],batch_visual_norm,tmp=TMP)
                all_cls_weights = batch_weights
                generalpreds2 = self.apply_classification_weights(batch_visual_norm, all_cls_weights)
                generalpreds1 = torch.argmax(generalpreds1,dim=1)
                generalpreds2 = torch.argmax(generalpreds2,dim=1)+len(self.seen_labels)
                generalpreds = generalpreds1
                for i in range(len(generalpreds1)):
                    if generalpreds1[i] > len(self.seen_labels)-1: # judge by model as unseen class
                        generalpreds[i] = generalpreds2[i]

                if dataType.lower()=='seen':
                    preds,_ = self.model1(self.attMatrix[:len(self.seen_labels)],images,TrainOrTest="Test")
                    preds = torch.argmax(preds,dim=1)
                elif dataType.lower()=='unseen':
                    preds = generalpreds2

                # seen or unseen split accuracy
                splitacc = [ i in seenOrunseenCls for i in generalpreds1.cpu().detach().tolist()]
                splitaccs.extend(splitacc)

                # calculate accs
                acc = preds == targets
                TF_l.extend(acc.cpu().detach().tolist())
                max_cos,_ = torch.max(generalpreds1_n,dim=-1)
                cos.extend(max_cos.cpu().detach().tolist())

                # calculate accs_cls
                for c in torch.unique(targets):
                    loc = targets == c
                    acc_c = preds[loc] == targets[loc]
                    generalacc_c = generalpreds[loc] == targets[loc]
                    dist_c = generalpreds1_n[loc,c].cpu().detach().tolist()
                    acc_c = acc_c.cpu().detach().tolist()
                    generalacc_c = generalacc_c.cpu().detach().tolist()
                    try:
                        accs_cls[c.cpu().detach().tolist()].extend(acc_c)
                        generalAccs_cls[c.cpu().detach().tolist()].extend(generalacc_c)
                        cos_cls[c.cpu().detach().tolist()].extend(dist_c)
                    except(KeyError):
                        accs_cls[c.cpu().detach().tolist()] = acc_c
                        generalAccs_cls[c.cpu().detach().tolist()] = generalacc_c
                        cos_cls[c.cpu().detach().tolist()] = dist_c

        # average class accuracy
        if accType == 'allCls':
            for c in accs_cls:
                accs.append(np.mean(accs_cls[c]))
                generalAccs.append(np.mean(generalAccs_cls[c]))
        for c in accs_cls:
            accs_cls[c] = np.mean(accs_cls[c])*100
            generalAccs_cls[c] = np.mean(generalAccs_cls[c])*100

        return  {'accs':np.mean(accs)*100,
                 'generalAccs':np.mean(generalAccs)*100,
                 'generalAccs_cls':generalAccs_cls,
                 'accs_cls':accs_cls,
                 'splitacc':np.mean(splitaccs)*100,
                 'cos_cls':cos_cls,
                 'TrueFalseList':TF_l,
                 'cos':cos
                }

    def prepareModel(self,):
        self.createLinearModel()

    def evalModel(self,chkFile1,chkFile2,resume=False):
        """
            train Model #epochs defined in cfg

            args:
                chkFile: save model chkpoint
        """

        testSeenData = self.loadDataFromCustom(dataType='test_seen')
        testUnseenData = self.loadDataFromCustom(dataType='test_unseen')

        self.loadChk(chkFile1,chkFile2)

        testSeenEval = self.evaluationSplit(testSeenData,dataType='seen',accType=self.cfg['model_hyp']['accType'],imgfType=self.cfg['dataset_hyp']['imgfType'],seenOrunseenCls=np.arange(0,len(self.seen_labels)))
        testUnseenEval = self.evaluationSplit(testUnseenData,dataType='unseen',accType=self.cfg['model_hyp']['accType'],imgfType=self.cfg['dataset_hyp']['imgfType'],seenOrunseenCls=np.arange(len(self.seen_labels),len(self.seen_labels)+len(self.unseen_labels)))

        accu = testUnseenEval['generalAccs']
        accs = testSeenEval['generalAccs']
        H = 2*(accu*accs)/(accu+accs)
        print('ZSL(T1):',testUnseenEval['accs'])   # ZSL(T1)
#         print('test/SeenAcc',testSeenEval['accs'])
        print('GZSL H:',H)                                # H
        print('GZSL Seen:',accs)                   # GZSL Seen
        print('GZSL Unseen:',accu)                   # GZSL Unseen
#         print('test/fu',testUnseenEval['splitacc'])
#         print('test/fs',testSeenEval['splitacc'])
        return testSeenEval['TrueFalseList'],testSeenEval['cos']

if __name__ =="__main__":

    parser = argparse.ArgumentParser()
    parser.add_argument('--modelChk',type=str,help='model chkpoint file')
    parser.add_argument('--yaml',type=str,help='checkpoint file saved in yaml')
    parser.add_argument('--device',default='cuda:0',type=str,help='cuda:0,cuda:1,cuda:2,cpu'
)
    opt = parser.parse_args()

    #chkFile1 ='./chk/201115CUBS2IFinetune_ResNet101_fullv_r520_IAStmp10CustomedExtractedCos520LR5e-4WD5e-7/chk_best_Hs.pt' #class weight + IAS
    chkFile1 ='./chk/CUBLR5e-4WD9e-7/chk_best_Hs.pt' #class weight + IAS
    chkFile2 = './chk/cublr5e-4_opt5e-7_w20_s4_51FinetunebestunseenAcc.pt'#mixup_IAS
    device = torch.device(opt.device)

    t = train(device,opt.yaml,IAS=True)
    t.prepareModel()

    truefalselistBase,cos_Base = t.evalModel(chkFile1,chkFile2)

reason: import numpy,from scipy,import networkx


text:
import pathlib
import os

import numpy as np

home = os.getenv('PT_DATA_DIR')
if home is None:
    home = str(pathlib.Path.home())
WEIGHTS_DIR = home + '/sm_weights/'

#  It's helpful to separate out the twenty canonical amino acids from the rest
CAN_AAS = 'ACDEFGHIKLMNPQRSTVWY'
AMB_AAS = 'BZX'
OTHER_AAS = 'JOU'
ALL_AAS = CAN_AAS + AMB_AAS + OTHER_AAS

DNA = 'GATC'
EXTENDED_NA = 'RYWSMKHBVDN'
RNA = 'GAUC'
IUPAC_AMB_DNA = DNA + EXTENDED_NA
IUPAC_AMB_RNA = RNA + EXTENDED_NA
NAS = 'GATUC' + EXTENDED_NA

STOP = '*'
GAP = '-'
PAD = GAP
MASK = '#'  # Useful for masked language model training
START = '@'

SPECIALS = STOP + GAP + MASK + START
PROTEIN_ALPHABET = ALL_AAS + SPECIALS
RNA_ALPHABET = IUPAC_AMB_RNA + SPECIALS

trR_ALPHABET = "ARNDCQEGHILKMFPSTWYV-"

AAINDEX_ALPHABET = 'ARNDCQEGHILKMFPSTWYV'

IUPAC_SS = 'HSTC'
DSSP = 'GHITEBSC'
SS8 = DSSP
SS3 = 'HSL'  # H: GHI; S: EB; L: STC

# Bins from TrRosetta paper
DIST_BINS = np.concatenate([np.array([np.nan]), np.linspace(2, 20, 37)])
THETA_BINS = np.concatenate([np.array([np.nan]), np.linspace(0, 2 * np.pi, 25)])
PHI_BINS = np.concatenate([np.array([np.nan]), np.linspace(0, np.pi, 13)])
OMEGA_BINS = np.concatenate([np.array([np.nan]), np.linspace(0, 2 * np.pi, 25)])

IUPAC_CODES = {
    "Ala": "A",
    "Arg": "R",
    "Asn": "N",
    "Asp": "D",
    "Cys": "C",
    "Gln": "Q",
    "Glu": "E",
    "Gly": "G",
    "His": "H",
    "Ile": "I",
    "Leu": "L",
    "Lys": "K",
    "Met": "M",
    "Phe": "F",
    "Pro": "P",
    "Ser": "S",
    "Thr": "T",
    "Trp": "W",
    "Val": "V",
    "Tyr": "Y",
    "Asx": "B",
    "Sec": "U",
    "Xaa": "X",
    "Glx": "Z",
}

max_stars_repo_name: microsoft/protein-sequence-models

max_issues_repo_name: microsoft/protein-sequence-models

max_forks_repo_name: microsoft/protein-sequence-models

reason: import numpy


text:
#encoding:UTF-8
#author:justry

import os
import numpy as np

convert_dict = {'0':'o',
                'o':'0',
                'O':'0',
                '1':'l',
                'l':'1'}

def ocr_correct(s):
    return s

def isalpha(letter):
    return letter.isalpha()

def isdigit(letter):
    return letter.isdigit()

def and_(a, b):
    return a and b

def or_(a, b):
    return a or b

def convert(i, s, func, l, boolean):
    if i == 0:
        if func(s[1]):
            l.append(i)
    elif i == len(s)-1:
        if func(s[-2]):
            l.append(i)
    elif boolean(func(s[i-1]), func(s[i+1])):
        l.append(i)

def o_0(s):
    pos = [i for i, j in enumerate(s) if j in '0oOl1']
    convert_list = []
    for i in pos:
        if s[i] in 'oOl':
            # letter o or O or l
            convert(i, s, isdigit, convert_list, and_)
        else:
            # digit 0, 1
            convert(i, s, isalpha, convert_list, or_)
    s = [i for i in s]
    for i in convert_list:
        s[i] = convert_dict[s[i]]
    return ''.join(s)

def has_end_char(s):
    return '.' in s or '!' in s or '?' in s or ';' in s

def uppercase_lowercase(s):
    s = [i for i in s.split(' ') if i != '']
    end_char_list = [has_end_char(i) for i in s]
    convert_list = []
    for i, j in enumerate(s):
        if i == 0:
            s[i] = s[i][0].upper() + s[i][1:].lower()
        else:
            if end_char_list[i-1]:
                s[i] = s[i][0].upper() + s[i][1:].lower()
            else:
                s[i] = s[i].lower()
    return ' '.join(s)

if __name__ == '__main__':
    file_path = '/home/zhaoke/justrypython/utils/data/test/9781601988157-summary14.result'
    with open(file_path, 'r') as f:
        lines = f.readlines()
    for i in lines:
        i = o_0(i)
        i = uppercase_lowercase(i)
        print(i)

max_stars_repo_name: justrypython/utils

max_issues_repo_name: justrypython/utils

max_forks_repo_name: justrypython/utils


text:
import matplotlib.pyplot as plt
import numpy as np

import pyk4a
from pyk4a import Config, PyK4A

MAX_SAMPLES = 1000

def set_default_data(data):
    data["acc_x"] = MAX_SAMPLES * [data["acc_x"][-1]]
    data["acc_y"] = MAX_SAMPLES * [data["acc_y"][-1]]
    data["acc_z"] = MAX_SAMPLES * [data["acc_z"][-1]]
    data["gyro_x"] = MAX_SAMPLES * [data["acc_x"][-1]]
    data["gyro_y"] = MAX_SAMPLES * [data["acc_y"][-1]]
    data["gyro_z"] = MAX_SAMPLES * [data["acc_z"][-1]]

def main():
    k4a = PyK4A(
        Config(
            color_resolution=pyk4a.ColorResolution.RES_720P,
            depth_mode=pyk4a.DepthMode.NFOV_UNBINNED,
            synchronized_images_only=True,
        )
    )
    k4a.start()

    plt.ion()
    fig, axes = plt.subplots(3, sharex=False)

    data = {
        "temperature": [0] * MAX_SAMPLES,
        "acc_x": [0] * MAX_SAMPLES,
        "acc_y": [0] * MAX_SAMPLES,
        "acc_z": [0] * MAX_SAMPLES,
        "acc_timestamp": [0] * MAX_SAMPLES,
        "gyro_x": [0] * MAX_SAMPLES,
        "gyro_y": [0] * MAX_SAMPLES,
        "gyro_z": [0] * MAX_SAMPLES,
        "gyro_timestamp": [0] * MAX_SAMPLES,
    }
    y = np.zeros(MAX_SAMPLES)
    lines = {
        "temperature": axes[0].plot(y, label="temperature")[0],
        "acc_x": axes[1].plot(y, label="acc_x")[0],
        "acc_y": axes[1].plot(y, label="acc_y")[0],
        "acc_z": axes[1].plot(y, label="acc_z")[0],
        "gyro_x": axes[2].plot(y, label="gyro_x")[0],
        "gyro_y": axes[2].plot(y, label="gyro_y")[0],
        "gyro_z": axes[2].plot(y, label="gyro_z")[0],
    }

    for i in range(MAX_SAMPLES):
        sample = k4a.get_imu_sample()
        sample["acc_x"], sample["acc_y"], sample["acc_z"] = sample.pop("acc_sample")
        sample["gyro_x"], sample["gyro_y"], sample["gyro_z"] = sample.pop("gyro_sample")
        for k, v in sample.items():
            data[k][i] = v
        if i == 0:
            set_default_data(data)

        for k in ("temperature", "acc_x", "acc_y", "acc_z", "gyro_x", "gyro_y", "gyro_z"):
            lines[k].set_data(range(MAX_SAMPLES), data[k])

        acc_y = data["acc_x"] + data["acc_y"] + data["acc_z"]
        gyro_y = data["gyro_x"] + data["gyro_y"] + data["gyro_z"]
        lines["acc_x"].axes.set_ylim(min(acc_y), max(acc_y))
        lines["gyro_x"].axes.set_ylim(min(gyro_y), max(gyro_y))
        lines["temperature"].axes.set_ylim(min(data["temperature"][0 : i + 1]), max(data["temperature"][0 : i + 1]))

        fig.canvas.draw()
        fig.canvas.flush_events()

    k4a._stop_imu()
    k4a.stop()

if __name__ == "__main__":
    main()


text:
from typing import List
import numpy as np
from pyrep.objects.shape import Shape
from pyrep.objects.proximity_sensor import ProximitySensor
from rlbench.const import colors
from rlbench.backend.task import Task
from rlbench.backend.conditions import DetectedCondition, ConditionSet
from rlbench.backend.spawn_boundary import SpawnBoundary

class HannoiSquare(Task):

    def init_task(self) -> None:
        self.square_ring = Shape('hannoi_square_ring')
        self.success_detector = ProximitySensor(
            'hannoi_square_success_detector')
        self.register_graspable_objects([self.square_ring])
        success_condition1 = DetectedCondition(self.square_ring,
                                               self.success_detector)
        self.register_success_conditions([success_condition1])

    def init_episode(self, index: int) -> List[str]:
        pillar0 = Shape('hannoi_square_pillar0')
        pillar1 = Shape('hannoi_square_pillar1')
        pillar2 = Shape('hannoi_square_pillar2')
        spokes = [pillar0, pillar1, pillar2]

        color_name, color_rgb = colors[index]
        chosen_pillar = np.random.choice(spokes)
        chosen_pillar.set_color(color_rgb)
        _, _, z = self.success_detector.get_position()
        x, y, _ = chosen_pillar.get_position()
        self.success_detector.set_position([x, y, z])

        color_choices = np.random.choice(
            list(range(index)) + list(range(index + 1, len(colors))),
            size=2, replace=False)
        spokes.remove(chosen_pillar)
        for spoke, i in zip(spokes, color_choices):
            name, rgb = colors[i]
            spoke.set_color(rgb)

        boundary1 = Shape('hannoi_square_boundary0')
        square_ring = Shape('hannoi_square_ring')
        b = SpawnBoundary([boundary1])
        b.sample(square_ring)

        return ['put the ring on the %s spoke' % color_name,
                'slide the ring onto the %s colored spoke' % color_name,
                'place the ring onto the %s spoke' %color_name]

    def variation_count(self) -> int:
        return len(colors)

max_stars_repo_name: TianhongDai/RLBench

max_issues_repo_name: TianhongDai/RLBench

max_forks_repo_name: TianhongDai/RLBench


text:
# -*- coding: utf-8 -*-
# -----------------------------------------------------------------------------
# MeshPy: A beam finite element input generator
#
# MIT License
#
# Copyright (c) 2021 Ivo Steinbrecher
#                    Institute for Mathematics and Computer-Based Simulation
#                    Universitaet der Bundeswehr Muenchen
#                    https://www.unibw.de/imcs-en
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.
# -----------------------------------------------------------------------------
"""
This script is used to test the functionality of the simulation manager.
"""

# Python imports.
import os
import unittest
import shutil
import numpy as np
import subprocess

# User imports.
from meshpy import (mpy, InputFile, MaterialReissner, BoundaryCondition,
    Beam3rHerm2Line3, Function, set_header_static, set_runtime_output)
from meshpy.mesh_creation_functions import create_beam_mesh_line
from meshpy.simulation_manager import (Simulation, SimulationManager)
from tests.testing_utility import (testing_temp, get_baci_path,
    testing_input, compare_strings)

def create_cantilever(convergence_base_dir, subpath, n_el):
    """
    Create a cantilever beam for a convergence analysis.
    """

    input_file = InputFile()
    set_header_static(input_file, time_step=0.25, n_steps=4)
    set_runtime_output(input_file, output_energy=True)
    mat = MaterialReissner(radius=0.1, youngs_modulus=10000.0)
    beam_set = create_beam_mesh_line(input_file, Beam3rHerm2Line3, mat,
        [0, 0, 0], [1, 0, 0], n_el=n_el)
    input_file.add(BoundaryCondition(beam_set['start'],
        ('NUMDOF 9 ONOFF 1 1 1 1 1 1 0 0 0 VAL 0 0 0 0 0 0 0 0 0 ' +
        'FUNCT 0 0 0 0 0 0 0 0 0'),
        bc_type=mpy.bc.dirichlet
        ))
    fun = Function('COMPONENT 0 FUNCTION t')
    input_file.add(fun, BoundaryCondition(beam_set['end'],
        ('NUMDOF 9 ONOFF 0 0 1 0 0 0 0 0 0 VAL 0 0 -{} 0 0 0 0 0 0 ' +
        'FUNCT 0 0 {} 0 0 0 0 0 0'),
        format_replacement=[0.5, fun],
        bc_type=mpy.bc.dirichlet
        ))
    base_dir = os.path.join(convergence_base_dir, subpath)
    if os.path.isdir(base_dir):
        shutil.rmtree(base_dir)
    os.makedirs(base_dir)
    input_file_path = os.path.join(base_dir, 'cantilever.dat')
    input_file.write_input_file(input_file_path, add_script_to_header=False,
        header=False)
    return input_file_path

class TestSimulationManager(unittest.TestCase):
    """Test various stuff of the simulation manager."""

    def setUp(self):
        """
        This method is called before each test and sets the default meshpy
        values for each test. The values can be changed in the individual
        tests.
        """

        # Set default values for global parameters.
        mpy.set_default_values()

    def create_simulations(self, convergence_base_dir):
        """
        Create the convergence study and return the SimulationManager
        """

        os.makedirs(convergence_base_dir, exist_ok=True)
        manager = SimulationManager(convergence_base_dir)
        manager.add(Simulation(
            create_cantilever(convergence_base_dir, 'ref', 40), n_proc=4))
        simulations = [Simulation(
            create_cantilever(os.path.join(convergence_base_dir, 'sim'),
                str(n_el), n_el)) for n_el
            in range(1, 7, 2)]
        manager.add(simulations)
        return manager

    def check_convergence_results(self, convergence_base_dir):
        """
        Check the results from the convergence study.
        """

        sim_dir = os.path.join(convergence_base_dir, 'sim')
        sub_dirs = [os.path.join(sim_dir, path) for path in os.listdir(sim_dir)
            if os.path.isdir(os.path.join(sim_dir, path))]
        sub_dirs.append(os.path.join(convergence_base_dir, 'ref'))
        results = {}
        for sub_dir in sub_dirs:
            my_data = np.genfromtxt(sub_dir + '/xxx_energy.csv',
                delimiter=',')
            key = sub_dir.split('/')[-1]
            results[key] = my_data[-1, 2]

        results_ref = {'5': 0.335081498526998, '3': 0.335055487040675,
            '1': 0.33453718896204, 'ref': 0.335085590674607}
        for key in results_ref.keys():
            self.assertTrue(abs(results[key] - results_ref[key]) < 1e-12)

    def xtest_simulation_manager(self, command):
        """
        Create a convergence study and check the results.

        Args
        ----
        command: str
            Name of the command that is used to run the study.
        """

        # Create the simulations.
        convergence_base_dir = os.path.join(testing_temp, 'simulation_manager')
        manager = self.create_simulations(convergence_base_dir)

        if shutil.which(command) is not None:
            if command == 'mpirun':
                manager.run_simulations_and_wait_for_finish(
                    baci_build_dir=os.path.dirname(get_baci_path()),
                    status=False)
            else:
                manager.submit_batch_files_and_wait_for_finish(
                    baci_build_dir=os.path.dirname(get_baci_path()),
                    check_interval=1)
        else:
            self.skipTest('{} was not found'.format(command))

        # Check the results.
        self.check_convergence_results(convergence_base_dir)

    def test_simulation_manager_mpirun(self):
        """
        Create a convergence study and check the results. The simulations are
        run with mpirun.
        """
        self.xtest_simulation_manager('mpirun')

    def test_simulation_manager_cluster(self):
        """
        Create a convergence study on the cluster and check the results.
        """
        self.xtest_simulation_manager('sbatch')

    def test_batch_file(self):
        """
        Test that the created batch file is correct.
        """

        sim = Simulation(os.path.join(testing_temp, 'dummy.dat'),
            n_proc=37,
            n_nodes=69,
            exclusive=True,
            output_prefix='xxxXXXxxx',
            wall_time='66:66:66',
            restart_step=17,
            restart_dir='../old_sim',
            restart_from_prefix='xxx_old',
            job_name='awsome_job',
            feature='skylake')

        sim.create_batch_file(testing_temp, 'batch_name')

        compare_strings(self,
            'test_batch_file',
            os.path.join(testing_temp, 'batch_name'),
            os.path.join(testing_input, 'batch_name_reference'))

if __name__ == '__main__':
    # Execution part of script.
    unittest.main()

max_stars_repo_name: imcs-compsim/meshpy

max_issues_repo_name: Manfred-Hyt/meshpy

max_forks_repo_name: Manfred-Hyt/meshpy


text:
from collections import OrderedDict

import numpy as np
import tensorflow as tf
from tensorflow.contrib.staging import StagingArea

from queue import deque
from baselines import logger
from baselines.template.util import logger as log_formater
from baselines.util import (
    import_function, store_args, flatten_grads, transitions_in_episode_batch)
from baselines.herhrl.normalizer import Normalizer
from baselines.herhrl.replay_buffer import ReplayBuffer
from baselines.common.mpi_adam import MpiAdam
from baselines.template.policy import Policy
from baselines.herhrl.hrl_policy import HRL_Policy
from baselines.herhrl.ddpg_her_hrl_policy import DDPG_HER_HRL_POLICY
from baselines.herhrl.pddl_policy import PDDL_POLICY
import math

def dims_to_shapes(input_dims):
    return {key: tuple([val]) if val > 0 else tuple() for key, val in input_dims.items()}

# class FlexiAvgList:
#     def __init__(self, maxlen=1000, n_chunks=4):
#         self.counter = 0
#         self.buffer = []
#         self.this_counter = 0
#         self.maxlen = maxlen
#
#     def add(self, item):
#         self.buffer.append(item)
#         self.counter += 1
#         self.this_counter += 1
#         if self.this_counter > self.maxlen:

class MIX_PDDL_HRL_POLICY(DDPG_HER_HRL_POLICY, PDDL_POLICY):
    # Putting DPG_HER_HRL_POLICY before PDDL_POLICY in inheritance assures that the functions in that class have precedence over those in PDDL_POLICY (see https://stackoverflow.com/questions/21657822/python-and-order-of-methods-in-multiple-inheritance)
    @store_args
    def __init__(self, input_dims, buffer_size, hidden, layers, network_class, polyak, batch_size,
                 Q_lr, pi_lr, norm_eps, norm_clip, max_u, action_l2, clip_obs, scope, T,
                 rollout_batch_size, subtract_goals, relative_goals, clip_pos_returns, clip_return,
                 sample_transitions, gamma, reuse=False, **kwargs):
        """Implementation of DDPG that is used in combination with Hindsight Experience Replay (HER).

        Args:
            input_dims (dict of ints): dimensions for the observation (o), the goal (g), and the
                actions (u)
            buffer_size (int): number of transitions that are stored in the replay buffer
            hidden (int): number of units in the hidden layers
            layers (int): number of hidden layers
            network_class (str): the network class that should be used (e.g. 'baselines.her.ActorCritic')
            polyak (float): coefficient for Polyak-averaging of the target network
            batch_size (int): batch size for training
            Q_lr (float): learning rate for the Q (critic) network
            pi_lr (float): learning rate for the pi (actor) network
            norm_eps (float): a small value used in the normalizer to avoid numerical instabilities
            norm_clip (float): normalized inputs are clipped to be in [-norm_clip, norm_clip]
            max_u (float): maximum action magnitude, i.e. actions are in [-max_u, max_u]
            action_l2 (float): coefficient for L2 penalty on the actions
            clip_obs (float): clip observations before normalization to be in [-clip_obs, clip_obs]
            scope (str): the scope used for the TensorFlow graph
            T (int): the time horizon for rollouts
            rollout_batch_size (int): number of parallel rollouts per DDPG agent
            subtract_goals (function): function that subtracts goals from each other
            relative_goals (boolean): whether or not relative goals should be fed into the network
            clip_pos_returns (boolean): whether or not positive returns should be clipped
            clip_return (float): clip returns to be in [-clip_return, clip_return]
            sample_transitions (function) function that samples from the replay buffer
            gamma (float): gamma used for Q learning updates
            reuse (boolean): whether or not the networks should be reused
        """
        DDPG_HER_HRL_POLICY.__init__(self, input_dims, buffer_size, hidden, layers, network_class, polyak, batch_size,
                 Q_lr, pi_lr, norm_eps, norm_clip, max_u, action_l2, clip_obs, scope, T,
                 rollout_batch_size, subtract_goals, relative_goals, clip_pos_returns, clip_return,
                 sample_transitions, gamma, reuse=reuse, **kwargs)

        self.count_pddl = 0.0
        self.count_hrl = 0.0
        # self.p_threshold = kwargs['p_threshold']
        self.p_threshold = {'train': 1.0, 'test': 1.0}
        self.p_steepness = kwargs['p_steepness']
        self.train_min_mix_entries = self.buffer_size // 8
        self.train_min_mix_entries = 2004
        # self.min_mix_entries = 20
        self.succ_rate_history = {'train': deque(maxlen=self.buffer_size // 2),
                                  'test': deque(maxlen=self.buffer_size // 2)}

    def get_actions(self, o, ag, g, noise_eps=0., random_eps=0., use_target_net=False,
                    compute_Q=False, exploit=True, success_rate=1.):

        def sigmoid(x):
            return 1 / (1 + math.exp(-x))

        def sigm_prob(x, p_threshold, steepness):
            # Verify using https://www.desmos.com/calculator 1/(1+\exp(-(2\cdot\left((x+0.5-0.3)\cdot12-6)))\right)
            offset_x = x - p_threshold
            scaled_x = offset_x * 12
            p = sigmoid(steepness * scaled_x)
            return p
        mode = 'train'
        if exploit:
            mode = 'test'

        self.succ_rate_history[mode].append(np.mean(success_rate))

        # If last 4/4 of history is not 5% better than last 3/4 of history after at least self.min_mix_entries entries, set switch point.
        # The switch point sets the threshold for sigmoidal policy sampling.
        if len(self.succ_rate_history['train']) > self.train_min_mix_entries and self.p_threshold[mode] == 1.0 and len(self.succ_rate_history['test']) > 4:
            frac_n = len(self.succ_rate_history[mode]) // 4
            prev = list(self.succ_rate_history[mode])[-(2*frac_n):-(frac_n)]
            last = list(self.succ_rate_history[mode])[-frac_n:]
            avg_prev = np.mean(prev)
            avg_last = np.mean(last)
            if avg_prev * 1.02 > avg_last and avg_last > 0.01:
                self.p_threshold[mode] = avg_last

        p_ddpg = sigm_prob(success_rate, self.p_threshold[mode], self.p_steepness)
        rnd = np.random.uniform()
        u, q = DDPG_HER_HRL_POLICY.get_actions(self, o, ag, g, noise_eps=noise_eps, random_eps=random_eps,
                                               use_target_net=use_target_net,
                                               compute_Q=compute_Q, exploit=exploit)
        if rnd > p_ddpg:
            scaled_u, q_pddl = PDDL_POLICY.get_actions(self, o, ag, g)
            # This u comes already in scale and with offsets, i.e., not in [-1,1]. So we have to descale it.
            u = self.inverse_scale_and_offset_action(scaled_u)
            if u.shape != g.shape:
                u = np.reshape(u, g.shape)  # this to solve the unbalance issue when rollout_batch_size = 1
            self.count_pddl += 1
        else:
            self.count_hrl += 1

        return u, q

    def logs(self, prefix='policy'):
        logs = []
        logs += [('ddpg_per_total', float(self.count_hrl / (self.count_hrl + self.count_pddl)))]
        logs += [('pddl_per_total', float(self.count_pddl / (self.count_hrl + self.count_pddl)))]
        logs += [('p_threshold_train', float(self.p_threshold['train']))]
        logs += [('p_threshold_test', float(self.p_threshold['test']))]

        logs = log_formater(logs, prefix + "_{}".format(self.h_level))
        logs += DDPG_HER_HRL_POLICY.logs(self, prefix=prefix)
        self.reset_counters()
        return logs

    def reset_counters(self):
        self.count_hrl = 0.0
        self.count_pddl = 0.0


text:
import numpy as np
import matplotlib.pyplot as plt  # To visualize
import pandas as pd  # To read data
from sklearn.linear_model import LinearRegression

data = pd.read_csv('fiyatlar.csv')  # veriyi alıyoruz
X = data.iloc[:, 0].values.reshape(-1, 1)  # Verideki değerleri numpy dizisine döndürüyoruz
Y = data.iloc[:, 2].values.reshape(-1, 1)  # -1, satırların boyutunun hesaplandığı, ancak 1 sütun olduğu anlamına gelir

linear_regressor = LinearRegression()  # model objemizi olusturduk
linear_regressor.fit(X, Y)  # lineer regresyon modelimizi eğittik
Y_pred = linear_regressor.predict(X)  # tahminler yaptık
print("----------------------------")
print(data)
print("--------------------------")
print(X)#eğittiğimiz x değerleri
print("--------------------------")
print(Y_pred)#tahminler
print("--------------------------")
print("Skor: " + str(linear_regressor.score(X,Y)))

plt.xlabel("Km ler")
plt.ylabel("Tahmini fiyatlar")
plt.scatter(X, Y)#x ve y değerlerine göre noktalarımızı grafiğe yerleştirdik.
plt.plot(X, Y_pred, color='red')
plt.show()

max_stars_repo_name: mertatay-ce/regression-analysis-project

max_issues_repo_name: mertatay-ce/regression-analysis-project

max_forks_repo_name: mertatay-ce/regression-analysis-project


text:
import torch
import torch.nn as nn
import torch.nn.functional as F
from torch.nn.utils.rnn import PackedSequence

import numpy as np

class SNN(nn.Module):
  @staticmethod
  def init(m):
    if type(m) == nn.Linear:
      nn.init.normal_(m.weight, std=np.sqrt(1./m.in_features))
      nn.init.constant_(m.bias, 0.)

  def __init__(self, layer_sizes):
    super(SNN, self).__init__()
    seq = []

    for i, d in enumerate(zip(layer_sizes[:-1], layer_sizes[1:])):
      if i > 0:
        seq.append(nn.SELU())
      seq.append(nn.Linear(*d))

    self.model = nn.Sequential(*seq)
    self.model.apply(SNN.init)

  def forward(self, input):
    return self.model(input)

class MelvinNet(nn.Module):
  def __init__(self, embedding_dim, hidden_dim, vocab_size, ff1_size, snn_hidden_dim, output_dim):
    super(MelvinNet, self).__init__()
    self.hidden_dim = hidden_dim
    self.embedding = nn.Embedding(vocab_size, embedding_dim)
    self.lstm = nn.LSTM(embedding_dim, hidden_dim)
    self.batchnorm = nn.BatchNorm1d(hidden_dim+ff1_size, affine=False)
    self.snn = SNN([hidden_dim+ff1_size] + snn_hidden_dim + [output_dim])

    # deactivate forget gate by initialization
    self.lstm.bias_hh_l0.data[hidden_dim:2*hidden_dim].fill_(float(np.finfo(np.float32).max))

  def init_hidden(self, batch_size, device):
    # The axes semantics are (num_layers, batch_size, hidden_dim)
    return (torch.zeros(1, batch_size, self.hidden_dim).to(device),
            torch.zeros(1, batch_size, self.hidden_dim).to(device))

  def forward(self, com, ff1, hidden):
    embeds = PackedSequence(self.embedding(com.data), com.batch_sizes)
    lstm_out, hidden = self.lstm(embeds, hidden)
    cat = torch.cat([hidden[1][0], ff1], -1)
    bn = self.batchnorm(cat)
    return self.snn(bn)

  def forward_lstm(self, com, hidden):
    embeds = PackedSequence(self.embedding(com.data), com.batch_sizes)
    lstm_out, hidden = self.lstm(embeds, hidden)
    return hidden[1][0]

  def forward_snn(self, ff1, cell):
    cat = torch.cat([cell, ff1], -1)
    bn = self.batchnorm(cat)
    return self.snn(bn)

class MelvinLstm(nn.Module):
  def __init__(self, vocab_size, embedding_dim, hidden_dim, output_dim):
    super(MelvinLstm, self).__init__()
    self.hidden_dim = hidden_dim
    self.embedding = nn.Embedding(vocab_size, embedding_dim)
    self.lstm = nn.LSTM(embedding_dim, hidden_dim)
    self.linear = nn.Linear(hidden_dim, output_dim)

    # init weights
    for i in range(4):
      torch.nn.init.orthogonal_(self.lstm.weight_hh_l0[hidden_dim * i : hidden_dim * (i+1)])
      torch.nn.init.xavier_uniform_(self.lstm.weight_ih_l0[hidden_dim * i : hidden_dim * (i+1)])

    # init biases to zero, deactivate forget gate
    self.lstm.bias_hh_l0.data.fill_(0.)
    self.lstm.bias_ih_l0.data.fill_(0.)
    self.lstm.bias_hh_l0.data[hidden_dim:2*hidden_dim].fill_(float(np.finfo(np.float32).max))
    self.lstm.bias_ih_l0.data[hidden_dim:2*hidden_dim].fill_(float(np.finfo(np.float32).max))

    self.linear.bias.data.fill_(0.)
    torch.nn.init.xavier_uniform_(self.linear.weight)

  def init_hidden(self, batch_size, device):
    # The axes semantics are (num_layers, batch_size, hidden_dim)
    return (torch.zeros(1, batch_size, self.hidden_dim).to(device),
            torch.zeros(1, batch_size, self.hidden_dim).to(device))

  def forward(self, com, hidden):
    embeds = PackedSequence(self.embedding(com.data), com.batch_sizes)
    lstm_out, hidden = self.lstm(embeds, hidden)
    return self.linear(hidden[0][0])

class Merlin(nn.Module):
  def __init__(self, io_size, hidden_size):
    super(Merlin, self).__init__()
    self.io_size = io_size
    self.hidden_size = hidden_size
    self.lstm = nn.LSTM(io_size, hidden_size)
    self.linear = nn.Linear(hidden_size, io_size)

    for i in range(4):
      torch.nn.init.orthogonal_(self.lstm.weight_hh_l0[hidden_size * i : hidden_size * (i+1)])
      torch.nn.init.xavier_uniform_(self.lstm.weight_ih_l0[hidden_size * i : hidden_size * (i+1)])

    self.lstm.bias_hh_l0.data.fill_(0.)
    self.lstm.bias_ih_l0.data.fill_(0.)
    self.lstm.bias_hh_l0.data[hidden_size:2*hidden_size].fill_(1.)
    self.lstm.bias_ih_l0.data[hidden_size:2*hidden_size].fill_(1.)
    torch.nn.init.xavier_uniform_(self.linear.weight)
    self.linear.bias.data.fill_(0.)

  def forward(self, com, hc):
    out, hc = self.lstm(com, hc)

    if isinstance(out, PackedSequence):
      out = PackedSequence(self.linear(out.data), out.batch_sizes)
    else:
      out = self.linear(out)

    return out, hc

max_stars_repo_name: ml-jku/melvin

max_issues_repo_name: ml-jku/melvin

max_forks_repo_name: ml-jku/melvin


text:
import policy as policy_module
import env_rna
import torch
import numpy as np
import torch.nn.functional as F
import matplotlib.pyplot as mlp
import arc_diagram
from torch.distributions import Categorical

env = env_rna.EnvRNA()

class Reinforce:

    running_reward = 0
    MAX_ITER = 1000
    exploration_eps = 10
    N = 10
    alpha = .90
    correct_predictions = 0
    sum_iterations_done = 0
    number_episodes = 100;
    policy = policy_module.Policy(N)
    optimizer = torch.optim.Adam(policy.parameters(), lr=1e-6)
    weights_dir = "./Weights/"
    def select_action(self,state, policy):
        probs = policy.forward(state)
        m = Categorical(probs)
        action = m.sample()
        policy.saved_probs.append(m.log_prob(action))
        item = action.item()
        return (int(item/self.N),item%self.N)

class MonteCarloReinforceTrainer(Reinforce):

    def train(self):

        for i_episode in range(self.number_episodes):
            if i_episode%10 == 0 :
                print("Episode : ", i_episode)
            seq = generate_random_sequence(self.N)
            env.reset(seq)
            state = env.rna.structure_representation
            ep_reward = 0
            rewards = []
            bestState = env.rna.structure_representation_dot.copy()
            bestReward = 0
            for t in range(self.MAX_ITER):
                exploration = np.random.uniform(0,1)
                if exploration > 0.3 or i_episode > self.exploration_eps or t == 0:#.9 * 1/(i_episode+1):
                    action = self.select_action(convert_to_tensor(state, seq), self.policy)
                    state, reward, done, _ = env.step(action,self.N)
                    rewards.append(reward)
                else:
                    action = np.random.randint(0,self.N,2)
                    action = (action[0],action[1])
                    state, reward, done, _ = env.step(action, self.N)
                    rewards.append(reward)
                ep_reward += reward
                if ep_reward > bestReward:
                    bestState = env.rna.structure_representation_dot.copy()
                    bestReward = ep_reward

                if done:
                    if i_episode > self.exploration_eps:
                        self.correct_predictions +=1
                        self.sum_iterations_done += t
                    print("Done ", i_episode, " iteration ", t)
                    title = str(i_episode) + " Done at iteration " + str(t)
                    if i_episode >= 2090:
                        mlp.show(arc_diagram.arc_diagram(
                            arc_diagram.phrantheses_to_pairing_list(env.rna.structure_representation_dot),seq, title))

                    break
                if (t+1)%100 == 0 and i_episode >= 2090:
                   mlp.show(arc_diagram.arc_diagram(arc_diagram.phrantheses_to_pairing_list(env.rna.structure_representation_dot),seq,i_episode))

            self.running_reward = self.running_reward * self.alpha + ep_reward * (1-self.alpha)
            if i_episode >= 2000:
                mlp.show(
                    arc_diagram.arc_diagram(arc_diagram.phrantheses_to_pairing_list(bestState),
                                            seq, "Best State Achieved"))

            self.finish_episode(rewards)
        self.policy.save_weights(self.weights_dir+ "monte_carlo_reinforce"+str(self.number_episodes))

    def finish_episode(self, rewards):
        R = 0
        DISCOUNT_FACTOR = 0.99
        policy_loss = []
        returns = []
        for r in rewards[::-1]:
            R = r+DISCOUNT_FACTOR*R
            returns.insert(0, R)
        returns = torch.tensor(returns)
        returns = (returns - returns.mean())
        for log_prob, R in zip(self.policy.saved_probs, returns):
            loss = -log_prob * R
            policy_loss.append(loss.unsqueeze(0))
        self.optimizer.zero_grad()
        policy_loss = torch.stack(policy_loss).sum()
        policy_loss.backward()
        self.optimizer.step()
        del self.policy.saved_probs[:]

class TemporalDifferenceReinforceTrainer(Reinforce):

    def train(self):

        for i_episode in range(self.number_episodes):
            seq = generate_random_sequence(self.N)
            env.reset(seq)
            state = env.rna.structure_representation
            ep_reward = 0
            rewards = []
            bestState = env.rna.structure_representation_dot.copy()
            bestReward = 0
            for t in range(self.MAX_ITER):
                exploration = np.random.uniform(0,1)
                if exploration > 0.3 or i_episode > 100 or t == 0:
                    action = self.select_action(convert_to_tensor(state, seq), self.policy)
                    state, reward, done, _ = env.step(action,self.N)
                    rewards.append(reward)
                    self.update_policy(reward)
                else:
                    action = np.random.randint(0,self.N,2)
                    action = (action[0],action[1])
                    state, reward, done, _ = env.step(action, self.N)
                    rewards.append(reward)

                ep_reward += reward
                if ep_reward > bestReward:
                    bestState = env.rna.structure_representation_dot.copy()
                    bestReward = ep_reward

                if done:
                    if i_episode > self.exploration_eps:
                        self.correct_predictions +=1
                        self.sum_iterations_done += t
                    print("Done ", i_episode, " iteration ", t)
                    title = str(i_episode) + " Done at iteration " + str(t)
                    if i_episode >= 3000:
                        mlp.show(arc_diagram.arc_diagram(
                            arc_diagram.phrantheses_to_pairing_list(env.rna.structure_representation_dot), seq, title))

                    break

                if (t + 1) % 100 == 0 and i_episode >= 3000:
                    mlp.show(arc_diagram.arc_diagram(
                        arc_diagram.phrantheses_to_pairing_list(env.rna.structure_representation_dot), seq, i_episode))

                self.running_reward = self.running_reward * self.alpha + ep_reward * (1 - self.alpha)
            if i_episode >= 3000:
                mlp.show(
                    arc_diagram.arc_diagram(arc_diagram.phrantheses_to_pairing_list(bestState),
                                            seq, "Best State Achieved"))

            self.finish_episode()
        self.policy.save_weights(self.weights_dir+"td_reinforce"+str(self.number_episodes))

    def update_policy(self,reward):
        R = 0
        DISCOUNT_FACTOR = 0.99
        policy_loss = []
        returns = []
        R = reward
        loss = -self.policy.saved_probs[-1] * R
        policy_loss.append(loss.unsqueeze(0))
        self.optimizer.zero_grad()
        policy_loss = F.smooth_l1_loss(2,loss)
        policy_loss.backward()
        self.optimizer.step()

    def finish_episode(self):
        del self.policy.saved_probs[:]

def convert_to_tensor(list,sequence):
    #TODO Change to NxN tensor
    n = len(sequence)
    tensor = torch.tensor(np.zeros((1,1,8,n)),dtype=torch.double)
    base_index = 0
    for base in sequence:
        position = 0
        if base == 'A': position = 0
        elif base == 'U': position = 2
        elif base == 'G' : position = 4
        else : position = 6
        if len([ (x,y) for x, y in list if x  == base_index or y == base_index ]) != 0:
            tensor[0][0][position+1][base_index] = 1
        else:
            tensor[0][0][position][base_index] = 1
        base_index +=1

    return tensor

def generate_random_sequence(N):
    sequence = ""
    for i in range (N):
        epsilon = np.random.uniform(0,1)
        if epsilon <= 0.25:
            sequence += "A"
        elif epsilon <= 0.5:
            sequence += "U"
        elif epsilon <= 0.75:
            sequence+= "G"
        else:
            sequence+= "C"

    return sequence

max_stars_repo_name: NetoPedro/RNAFoldingDeepRL

max_issues_repo_name: NetoPedro/RNAFoldingDeepRL

max_forks_repo_name: NetoPedro/RNAFoldingDeepRL


text:
from math import sqrt, ceil
import numpy as np
from pychemia.utils.serializer import PyChemiaJsonable, generic_serializer

"""
Definition of the set of k-points in reciprocal space
"""

__author__ = "Guillermo Avendano-Franco"
__copyright__ = "Copyright 2014"
__version__ = "1.1"
__maintainer__ = "Guillermo Avendano-Franco"
__email__ = "gtux.gaf@gmail.com"
__status__ = "Development"
__date__ = "November 5, 2014"

class KPoints(PyChemiaJsonable):
    """
    Defines an object that contains information about kpoints mesh
    There are two mutually exclusive options:

    Explicit array of kpoints with associated weights OR
    a specification of a grid with an optional set of shifts.

    The default is one single kpoint in gamma
    """

    def __init__(self, kmode='gamma',
                 kpoints_list=None, weights=None,
                 shifts=None, grid=None,
                 kvertices=None, intermediates=None):
        """
        The KPoints class is designed to  store the description of k-points in a code agnostic way.

        There are 5 possible kmodes to describe a set of k-points, currently those modes are not interchangable,
        but that could change in the future.

        The kmodes are 'cartesian', 'reduced', 'gamma', 'monkhorst-pack' and 'path'.

        Kmode 'cartesian' and 'reduced' uses explicit declaration of a grid of k-points.
        In the case of 'cartesian' the points are vectors in k-space.
        In the case of 'reduced' the points are the factors for the vectors defining the reciprocal lattice.
        The variable 'kpoints_list' is a numpy array with a shape of ('nkpt', 3)
        The variable 'weigths' is a numpy array with the weights associated to each point.
        by default all the weights are equal to  1.

        Kmode 'gamma' and 'monkhorst-pack' uses an uniform grid to describe the k-points mesh.
        The variable 'grid' describes the number of points on each dimension an it is a numpy array
        with 3 values.
        The variable 'shifts' describe the shift done to the grid according to the mode selected.
        The internal value 'nkpt' will be simply the product of the three values of the grid.
        Typically ab-initio packages will try to exclude from that grid equivalent k-points due to
        symmetry considerations.

        Kmode 'path' is specially useful for band-structure calculations, in this case the variable
        'kvertices' contains all the vertices, usually high symmetry points and the variable
        'intermediates' contains the number of k-points created between two consecutive vertices.

        When a Kpoints object is created without any arguments, the default will be a kmode='gamma',
        with a grid=[1, 1, 1] and shifts=[0, 0, 0]

        :param kmode: The options are 'cartesian', 'reduced', 'gamma', 'monkhorst-pack'
        :param kpoints_list: Explicit list of k-points (Only for 'cartesian' and 'reduced')
        :param weights: List of weights associated to each k-point (Only for 'cartesian' and 'reduced')
        :param grid: Number of kpoints on each direction (Only for 'gamma' and 'monkhorst-pack')
        :param shifts: Shift applied to the grid (Only for 'gamma' and 'monkhorst-pack')
        :param kvertices: List of vertices for a path of k-points (Only for 'path')
        :param intermediates: Number of intermediates between two vertices (Only for 'path')

        :rtype : Kpoints object
        """
        if kmode.lower() not in ['cartesian', 'reduced', 'gamma', 'monkhorst-pack', 'path']:
            raise ValueError("kmode must be 'cartesian', 'reduced', 'gamma', 'monkhorst-pack' or 'path'")

        self.kmode = kmode.lower()

        # All the variables initially None
        self.kpoints_list = None
        self.weights = None
        self.grid = None
        self.shifts = None
        self.kvertices = None
        self.intermediates = None

        # Explicit list of k-points
        if self.kmode in ['cartesian', 'reduced']:
            self.set_kpoints_list(kpoints_list, weights)
        # Grid specification
        elif self.kmode in ['gamma', 'monkhorst-pack']:
            self.set_grid(grid, shifts)
        # Path along the K-space
        elif self.kmode == 'path':
            if kvertices is None:
                self.kvertices = np.array([[0, 0, 0]])
            else:
                if np.array(kvertices).shape == (3,):
                    self.kvertices = np.array(kvertices).reshape((-1, 3))
                elif np.array(kvertices).shape[1] == 3:
                    self.kvertices = np.array(kvertices)
                else:
                    raise ValueError("Wrong value for kvertices, should be an array with shape (nkpt,3)")
            if intermediates is None:
                self.intermediates = 1
            else:
                self.intermediates = int(intermediates)

    def add_kpt(self, kpoint, weight=1):
        """
        Add one extra kpoint to the kpoints_list in the case of
        kmode 'cartesian' or 'reduced'

        :param kpoint: (list, numpy.array) The kpoint to be added
        :param weight: (int, float) the weight associated to that point
        """
        kpoint = list(kpoint)
        assert len(kpoint) == 3
        assert (self.kmode in ['cartesian', 'reduced'])
        if kpoint not in self.kpoints_list:
            self.kpoints_list.append(kpoint)
            self.weights = np.append(self.weights, weight)
        else:
            position = self.kpoints_list.index(kpoint)
            self.weights[position] = weight

    def set_kpoints_list(self, kpoints_list, weights=None):
        """
        Set an explicit list of kpoints with a proper check of correct arguments

        :param kpoints_list: Explicit list of k-points (Only for 'cartesian' and 'reduced')
        :param weights: List of weights associated to each k-point (Only for 'cartesian' and 'reduced')
        """
        assert (self.kmode in ['cartesian', 'reduced'])
        if kpoints_list is None:
            self.kpoints_list = [[0, 0, 0]]
        else:
            if np.array(kpoints_list).shape == (3,):
                self.kpoints_list = np.array(kpoints_list).reshape((-1, 3))
            elif np.array(kpoints_list).shape[1] == 3:
                self.kpoints_list = np.array(kpoints_list)
            else:
                raise ValueError("Wrong value for kpoints_list, should be an array with shape (nkpt,3)")
        self.kpoints_list = generic_serializer(self.kpoints_list)
        nkpt = len(self.kpoints_list)
        if weights is None:
            self.weights = np.ones(nkpt)
        elif len(np.array(weights)) == nkpt:
            self.weights = np.array(weights)
        else:
            raise ValueError("Wrong value for weights, should be an array with shape (nkpt,)")
        self.weights = list(self.weights)

    def set_grid(self, grid, shifts=None):
        """
        Set a grid of kpoints with a proper check of correct arguments

        :param grid: Number of kpoints on each direction (Only for 'gamma' and 'monkhorst-pack')
        :param shifts: Shift applied to the grid (Only for 'gamma' and 'monkhorst-pack')
        """

        assert (self.kmode in ['gamma', 'monkhorst-pack'])
        if grid is None:
            self.grid = [1, 1, 1]
        elif len(grid) == 3:
            self.grid = [int(x) for x in grid]
        else:
            raise ValueError("Wrong value for grid, should be an array with shape (-1,3)")
        if shifts is None:
            self.shifts = [0, 0, 0]
        else:
            self.shifts = generic_serializer(np.array(shifts).reshape(3))

    @property
    def number_of_kpoints(self):
        return self.nkpt

    @property
    def nkpt(self):
        """
        Number of kpoints
        For a grid the value is the product of the number of kpoints
        on each direction.
        """
        if self.kmode in ['cartesian', 'reduced']:
            return len(self.kpoints_list)
        elif self.kmode in ['gamma', 'monkhorst-pack']:
            return np.prod(self.grid)
        elif self.kmode == 'path':
            return len(self.kvertices) * self.intermediates

    def __str__(self):
        """
        String representation of the KPoints object
        """
        kp = ' Mode  : ' + self.kmode + '\n'
        if self.kmode in ['cartesian', 'reduced']:
            kp += ' Number of k-points: ' + str(self.nkpt) + '\n\n'
            for i in range(self.nkpt):
                kp += (" %10.5f %10.5f %10.5f %15.5f\n"
                       % (self.kpoints_list[i][0], self.kpoints_list[i][1], self.kpoints_list[i][2], self.weights[i]))
        elif self.kmode in ['gamma', 'monkhorst-pack']:
            kp += ' Grid  : ' + str(self.grid) + '\n'
            kp += ' Shift : ' + str(self.shifts) + '\n'
        elif self.kmode == 'path':
            kp += ' Number of intermediates: ' + str(self.intermediates) + '\n\n'
            for vertex in self.kvertices:
                kp += " %10.5f %10.5f %10.5f\n" % (vertex[0], vertex[1], vertex[2])
        return kp

    @property
    def to_dict(self):

        ret = {}
        for i in list(self.__dict__.keys()):
            if eval('self.' + i) is not None:
                ret[i] = eval('self.' + i)
        return ret

    @classmethod
    def from_dict(cls, json_dict):
        if json_dict['kmode'] in ['cartesian', 'reduced']:
            return cls(kmode=json_dict['kmode'], kpoints_list=json_dict['kpoints_list'],
                       weights=json_dict['weights'])
        elif json_dict['kmode'] in ['gamma', 'monkhorst-pack']:
            return cls(kmode=json_dict['kmode'], grid=json_dict['grid'], shifts=json_dict['shifts'])
        elif json_dict['kmode'] == 'path':
            return cls(kmode=json_dict['kmode'], kvertices=json_dict['kvertices'],
                       intermediates=json_dict['intermediates'])

    @staticmethod
    def max_kpoints(structure, kpoints_per_atom=1000):
        """
        Returns the maximal grid of kpoints allowed during the
        kpoints convergence test. Based on the routine used in
        pymatgen.

        :param structure: (pychemia.Structure) The structure to compute the maximal grid
        :param kpoints_per_atom: (int) number of kpoints per atom
        :return:
        """

        lengths = [sqrt(sum(map(lambda y: y ** 2, structure.cell[i]))) for i in range(3)]
        ngrid = kpoints_per_atom / structure.natom
        mult = (ngrid * lengths[0] * lengths[1] * lengths[2]) ** (1.0 / 3.0)
        num_div = [int(ceil(1.0 / lengths[i] * mult)) for i in range(3)]
        num_div = [i if i > 0 else 1 for i in num_div]
        return num_div[0], num_div[1], num_div[2]

    def get_density_of_kpoints(self, lattice):

        assert (self.kmode in ['gamma', 'monkhorst-pack'])
        vol = 1.0 / lattice.volume
        return self.number_of_kpoints / vol

    @classmethod
    def optimized_grid(cls, lattice, kp_density=1E5, nkpoints=None, force_odd=False):
        """
        Returns a KPoints object with kmode='gamma' with a grid optimized for a given
         density or number of kpoints.
         If force_odd is True the grid will be such that there are odd numbers on each direction
         such that 'gamma' is always included on the grid.

        :param lattice: (pychemia.Lattice) A PyChemia Lattice object, if the argument is a pychemia.Structure
                        object, its lattice will be used instead.
        :param kp_density: (float) A density of kpoints (default: 1E5)
        :param nkpoints: (int) Aproximate number of kpoints to use, if this argument is present the kp_density
                                will be ignored
        :param force_odd: (bool) If True, force the grid to be odd for the three directions.
        """
        if nkpoints is not None:
            kp_density = None

        # Reciprocal lattice object
        rlattice = lattice.reciprocal()

        if nkpoints is not None:
            vol = 1.0 / lattice.volume
            kp_density = nkpoints / vol

        factor = int(kp_density ** 1.0 / 3.0)

        while True:
            ini_density = kp_density

            grid = np.array([int(max(ceil(factor * rlattice.a), 1)),
                             int(max(ceil(factor * rlattice.b), 1)),
                             int(max(ceil(factor * rlattice.c), 1))])

            # The volume in reciprocal space is the inverse of the volume in real space
            fin_density = float(np.prod(grid) * lattice.volume)

            if fin_density > ini_density and factor > 0:
                factor -= 1
            else:
                break

        if force_odd:
            for i in range(3):
                if grid[i] % 2 == 0:
                    grid[i] += 1
        return cls(kmode='gamma', grid=grid)

max_stars_repo_name: petavazohi/PyChemia

max_issues_repo_name: petavazohi/PyChemia

max_forks_repo_name: petavazohi/PyChemia


text:
import numpy as np
import pyemma
from pyemma.util.contexts import settings
import MDAnalysis as mda
import matplotlib.pyplot as plt
from pensa.preprocessing import sort_coordinates, merge_and_sort_coordinates
from .visualization import project_on_eigenvector, sort_traj_along_projection

# --- METHODS FOR TIME-LAGGED INDEPENDENT COMPONENT ANALYSIS ---

# http://emma-project.org/latest/api/generated/pyemma.coordinates.tica.html

def calculate_tica(data, dim=-1, lag=10):
    """
    Performs a PyEMMA TICA on the provided data.

    Parameters
    ----------
        data : float array
            Trajectory data. Format: [frames,frame_data].
        dim : int, optional, default -1
            The number of dimensions (independent components) to project onto.
            -1 means all numerically available dimensions will be used.
        lag : int, optional, default = 10
            The lag time, in multiples of the input time step.

    Returns
    -------
        tica : TICA obj
            Time-lagged independent component information.

    """
    tica = pyemma.coordinates.tica(data, lag=lag)
    return tica

def tica_eigenvalues_plot(tica, num=12, plot_file=None):
    """
    Plots the highest eigenvalues over the number of the time-lagged independent components.

    Parameters
    ----------
        tica : TICA obj
            Time-lagged independent components information.
        num : int, default = 12
            Number of eigenvalues to plot.
        plot_file : str, optional, default = None
            Path and name of the file to save the plot.

    """
    # Plot eigenvalues over component numbers.
    fig,ax = plt.subplots(1, 1, figsize=[4,3], dpi=300)
    componentnr = np.arange(num)+1
    eigenvalues = tica.eigenvalues[:num]
    ax.bar(componentnr, eigenvalues)
    ax.set_xlabel('component number')
    ax.set_ylabel('eigenvalue')
    fig.tight_layout()
    # Save the figure to a file.
    if plot_file: fig.savefig(plot_file, dpi=300)
    return componentnr, eigenvalues

def tica_features(tica, features, num, threshold, plot_file=None, add_labels=False):
    """
    Prints relevant features and plots feature correlations.

    Parameters
    ----------
        tica : TICA obj
            The TICA of which to plot the features.
        features : list of str
            Features for which the TICA was performed
            (obtained from features object via .describe()).
        num : float
            Number of feature correlations to plot.
        threshold : float
            Features with a correlation above this will be printed.
        plot_file : str, optional, default = None
            Path and name of the file to save the plot.
        add_labels : bool, optional, default = False
            Add labels of the features to the x axis.

    """
    # Plot the highest TIC correlations and print relevant features.
    height = num*2+2 if add_labels else num*2
    fig,ax = plt.subplots(num,1,figsize=[4,height],dpi=300,sharex=True)
    for i in range(num):
        relevant = tica.feature_TIC_correlation[:,i]**2 > threshold**2
        print("Features with abs. corr. above a threshold of %3.1f for TIC %i:"%(threshold, i+1))
        for j, ft in enumerate(features):
            if relevant[j]: print(ft, "%6.3f"%(tica.feature_TIC_correlation[j,i]))
        ax[i].plot(tica.feature_TIC_correlation[:,i])
        test_feature = tica.feature_TIC_correlation[:,i]
        ax[i].set_ylabel('corr. with TIC%i'%(i+1))
    if add_labels:
        ax[-1].set_xticks(np.arange(len(features)))
        ax[-1].set_xticklabels(features,rotation=90)
    else:
        ax[-1].set_xlabel('feature index')
    fig.tight_layout()
    # Save the figure to a file.
    if plot_file: fig.savefig(plot_file,dpi=300)
    return test_feature

def project_on_tic(data, ev_idx, tica=None, dim=-1, lag=10):
    """
    Projects a trajectory onto an eigenvector of its TICA.

    Parameters
    ----------
        data : float array
            Trajectory data [frames,frame_data].
        ev_idx : int
            Index of the eigenvector to project on (starts with zero).
        tica : TICA obj, optional, default = None
            Information of pre-calculated TICA.
            Must be calculated for the same features (but not necessarily the same trajectory).
        dim : int, optional, default = -1
            The number of dimensions (independent components) to project onto.
            Only used if tica is not provided.
        lag : int, optional, default = 10
            The lag time, in multiples of the input time step.
            Only used if tica is not provided.

    Returns
    -------
        projection : float array
            Value along the TIC for each frame.

    """
    # Perform standard TICA if none is provided.
    if tica is None: calculate_tica(data, dim=dim, lag=lag)
    # Project the features onto the time-lagged independent components.
    projection = project_on_eigenvector(data, ev_idx, tica)
    return projection

def get_components_tica(data, num, tica=None, dim=-1, lag=10, prefix=''):
    """
    Projects a trajectory onto the first num eigenvectors of its tICA.

    Parameters
    ----------
        data : float array
            Trajectory data [frames,frame_data].
        num : int
            Number of eigenvectors to project on.
        tica : tICA obj, optional, default = None
            Information of pre-calculated tICA. Defaults to None.
            Must be calculated for the same features (but not necessarily the same trajectory).
        dim : int, optional, default = -1
            The number of dimensions (independent components) to project onto.
            Only used if tica is not provided.
        lag : int, optional, default = 10
            The lag time, in multiples of the input time step.
            Only used if tica is not provided.
        prefix : str, optional, default = ''
            First part of the component names. Second part is "IC"+<IC number>

    Returns
    -------
        comp_names : list
            Names/numbers of the components.
        components : float array
            Component data [frames,components]

    """
    # Perform tICA if none is provided.
    if tica is None: calculate_tica(data, lag=lag)
    # Project the features onto the principal components.
    comp_names = []
    components = []
    for ev_idx in range(num):
        projection = np.zeros(data.shape[0])
        for ti in range(data.shape[0]):
            projection[ti] = np.dot(data[ti],tica.eigenvectors[:,ev_idx])
        components.append(projection)
        comp_names.append(prefix+'IC'+str(ev_idx+1))
    # Return the names and data.
    return comp_names, np.array(components).T

def sort_traj_along_tic(data, top, trj, out_name, tica=None, num_ic=3, lag=10, start_frame=0):
    """
    Sort a trajectory along independent components.

    Parameters
    ----------
        data : float array
            Trajectory data [frames,frame_data].
        top : str
            File name of the reference topology for the trajectory.
        trj : str
            File name of the trajetory from which the frames are picked.
            Should be the same as data was from.
        out_name : str
            Core part of the name of the output files
        tica : tICA obj, optional, default = None
            Time-lagged independent components information.
            If none is provided, it will be calculated.
            Defaults to None.
        num_ic : int, optional, default = 3
            Sort along the first num_ic independent components.
        lag : int, optional, default = 10
            The lag time, in multiples of the input time step.
            Only used if tica is not provided.
        start_frame : int, optional, default = 0
            Offset of the data with respect to the trajectories (defined below).

    Returns
    -------
        sorted_proj: list
            sorted projections on each principal component
        sorted_indices_data : list
            Sorted indices of the data array for each principal component
        sorted_indices_traj : list
            Sorted indices of the coordinate frames for each principal component

    """
    # Calculate the principal components if they are not given.
    if tica is None: tica = calculate_tica(data, dim=num_ic, lag=lag)
    # Sort the trajectory along them.
    sorted_proj, sorted_indices_data, sorted_indices_traj = sort_traj_along_projection(
        data, tica, top, trj, out_name, num_comp=num_ic, start_frame=start_frame
    )
    return sorted_proj, sorted_indices_data, sorted_indices_traj

def sort_trajs_along_common_tic(data_a, data_b, top_a, top_b, trj_a, trj_b, out_name, num_ic=3, lag=10, start_frame=0):
    """
    Sort two trajectories along their most important common time-lagged independent components.

    Parameters
    ----------
        data_a : float array
            Trajectory data [frames,frame_data].
        data_b : float array
            Trajectory data [frames,frame_data].
        top_a : str
            Reference topology for the first trajectory.
        top_b : str
            Reference topology for the second trajectory.
        trj_a : str
            First of the trajetories from which the frames are picked.
            Should be the same as data_a was from.
        trj_b : str
            Second of the trajetories from which the frames are picked.
            Should be the same as data_b was from.
        out_name : str
            Core part of the name of the output files.
        num_ic : int, optional, default = 3
            Sort along the first num_ic independent components.
        lag : int, optional, default = 10
            The lag time, in multiples of the input time step.
            Only used if tica is not provided.
        start_frame : int, optional, default = 0
            Offset of the data with respect to the trajectories (defined below).

    Returns
    -------
        sorted_proj: list
            sorted projections on each principal component
        sorted_indices_data : list
            Sorted indices of the data array for each principal component
        sorted_indices_traj : list
            Sorted indices of the coordinate frames for each principal component

    """
    sorted_proj, sorted_indices_data, sorted_indices_traj = sort_mult_trajs_along_common_tic(
        [data_a, data_b], [top_a, top_b], [trj_a, trj_b], out_name, num_ic=num_ic, lag=lag, start_frame = start_frame
        )
    return sorted_proj, sorted_indices_data, sorted_indices_traj

def sort_mult_trajs_along_common_tic(data, top, trj, out_name, num_ic=3, lag=10, start_frame=0):
    """
    Sort multiple trajectories along their most important independent components.

    Parameters
    ----------
        data : list of float arrays
            List of trajectory data arrays, each [frames,frame_data].
        top : list of str
            Reference topology files.
        trj : list of str
            Trajetories from which the frames are picked.
            trj[i] should be the same as data[i] was from.
        out_name : str
            Core part of the name of the output files.
        num_ic : int, optional, default = 3
            Sort along the first num_ic independent components.
        lag : int, optional, default = 10
            The lag time, in multiples of the input time step.
            Only used if tica is not provided.
        start_frame : int or list of int, default = 0
            Offset of the data with respect to the trajectories.

    Returns
    -------
        sorted_proj: list
            sorted projections on each independent component
        sorted_indices_data : list
            Sorted indices of the data array for each independent component
        sorted_indices_traj : list
            Sorted indices of the coordinate frames for each independent component

    """
    num_frames = [len(d) for d in data]
    num_traj = len(data)
    if type(start_frame) == int:
        start_frame *= np.ones(num_traj)
        start_frame = start_frame.tolist()
    # Combine the input data
    all_data = np.concatenate(data,0)
    # Calculate the independent components
    tica = pyemma.coordinates.tica(all_data, lag=lag)
    # Initialize output
    sorted_proj = []
    sorted_indices_data = []
    sorted_indices_traj = []
    # Loop over principal components.
    for evi in range(num_ic):
        # Project the combined data on the independent component
        proj = [project_on_tic(d, evi, tica=tica) for d in data]
        # Sort everything along the projection on the respective independent component
        out_xtc = out_name+"_tic"+str(evi+1)+".xtc"
        proj_sort, sort_idx, oidx_sort = merge_and_sort_coordinates(
            proj, top, trj, out_xtc, start_frame=start_frame, verbose=False
            )
        sorted_proj.append(proj_sort)
        sorted_indices_data.append(sort_idx)
        sorted_indices_traj.append(oidx_sort)
    return sorted_proj, sorted_indices_data, sorted_indices_traj

max_stars_repo_name: drorlab/pensa

max_issues_repo_name: drorlab/pensa

max_forks_repo_name: drorlab/pensa


text:
# -*- coding: utf-8 -*-
"""
Created on Sat Apr  4 19:49:23 2020

@author: wmont
"""

import numpy as np

def calculate_fitness(element, matrix):
    '''Calculates the fitness of a given element based on the distances informed in the matrix.'''
    distance = 0
    for i in range(len(element)-1):
        distance += matrix[element[i]][element[i+1]]
    return distance

def initialize_population(n_pop, n_elements):
    '''Initializes a random population.'''
    population = [0]*n_pop

    list_elements = list(range(n_elements))
    for i in range(n_pop):
        # generating the element order
        element = list_elements.copy()
        np.random.shuffle(element)
        population[i] = element

    return population

def random_selection(population, fitness):
    '''Randomly chooses elements out of the population.'''
    normalized_fitness = fitness / np.sum(fitness)

    ranking = np.argsort(normalized_fitness)
    normalized_fitness = np.sort(normalized_fitness)
    cumulative_sum = np.cumsum(normalized_fitness)

    # since we want to minimize the distances, we're using argmin
    chosen_element = np.argmin(cumulative_sum >= np.random.random())
    return population[ranking[chosen_element]]

def generate_cuts(elements):
    '''Generates the cuts for the crossover.'''
    cuts = [0, 0]
    cuts[0] = np.random.randint(elements-2)+2

    while True:
        cuts[1] = np.random.randint(elements-2)+2
        if cuts[1] != cuts[0]:
            break

    return cuts

def remove_elements(elements, elements_to_ignore):
    '''Eliminates items from the first list contained in the second list.'''
    return [element for element in elements if element not in
            elements_to_ignore]

def crossover(population, new_population, crossover_size, crossover_rate):
    '''Crossover step.'''
    # based on lecture IS 07, slide 33 (order crossover)
    offspring = []
    n_elements = len(population[0])

    cuts = np.random.randint(n_elements-crossover_size-1)
    cuts = [cuts, cuts+5]

    for i in range(len(population)):
        # generating the offspring based on their parents
        offspring_1 = population[i].copy()
        offspring_2 = new_population[i].copy()

        if (np.random.random() <= crossover_rate):
            # copying the bits between the cutting points
            cut_1 = population[i][cuts[0]:cuts[1]]
            cut_2 = new_population[i][cuts[0]:cuts[1]]

            # getting the bit sequence from the second cut onwards
            order_1 = remove_elements(new_population[i][
                    cuts[1]:] + new_population[i][:cuts[1]], cut_1)

            offspring_1[cuts[1]:] = order_1[:n_elements-cuts[1]]
            offspring_1[:cuts[0]] = order_1[n_elements-cuts[1]:]

            order_2 = remove_elements(population[i][
                    cuts[1]:] + population[i][:cuts[1]], cut_2)
            offspring_2[cuts[1]:] = order_2[:n_elements-cuts[1]]
            offspring_2[:cuts[0]] = order_2[n_elements-cuts[1]:]

        offspring.append(offspring_1)
        offspring.append(offspring_2)

    return offspring

def mutation(population, mutation_rate):
    '''Mutation step.'''
    # based on lecture IS 07, slide 40 (reciprocal exchange)
    offspring = []
    n_elements = len(population[0])

    for i in range(len(population)):
        child = population[i].copy()
        if (np.random.random() <= mutation_rate):
            first_element = np.random.randint(n_elements-2)
            second_element = np.random.randint(first_element, n_elements-1)

            child[first_element] = population[i][second_element]
            child[second_element] = population[i][first_element]

        offspring.append(child)

    return offspring

def run_ga(matrix, pop, n_pop, max_iter, crossover_size, crossover_rate,
           mutation_rate, debug):
    '''Genetic Algorithm for Traveling Salesman Problems (TSP).'''
    # current iteration
    current_iteration = 1

    # while the stopping criteria had not been met
    while current_iteration <= max_iter:
        # fitness calculation
        fitness = [0] * n_pop
        for i in range(len(pop)):
            fitness[i] = calculate_fitness(pop[i], matrix)

        new_pop = pop.copy()
        for i in range(len(pop)):
            # selection
            new_pop[i] = random_selection(pop, fitness)

        # crossover
        cross_pop = crossover(pop, new_pop, crossover_size, crossover_rate)

        # mutation
        mut_pop = mutation(cross_pop, mutation_rate)

        # fitness calculation
        fitness = [0] * len(mut_pop)
        for i in range(len(mut_pop)):
            fitness[i] = calculate_fitness(mut_pop[i], matrix)

        # pruning (limiting the number of elements by n_pop)
        ranking = np.argsort(fitness)[:n_pop]
        fitness = list(fitness[i] for i in ranking)
        pop = list(mut_pop[i] for i in ranking)

        # showing stats
        if debug:
            print(f'Iteration {current_iteration}, ' +
                  f'min. distance {np.min(fitness)}, ' +
                  f'max. distance {np.max(fitness)}.')

        # incrementing the counter
        current_iteration += 1

    return fitness, pop