text:
#!/usr/bin/env python
#------------------------------------------------------------
# Script compares efficiency of automatic derivatives vs
# analytical in mpfit.py
# Vog, 31 okt 2011
#------------------------------------------------------------

import numpy
from matplotlib.pyplot import figure, show, rc
from kapteyn import kmpfit
from matplotlib.patches import Polygon

def confpred_band(x, dfdp, prob, fitobj, f, prediction, abswei=False, err=None):
   #----------------------------------------------------------
   # Return values for a confidence or a prediction band.
   # See documentation for methods confidence_band and
   # prediction_band
   #----------------------------------------------------------
   from scipy.stats import t
   # Given the confidence or prediction probability prob = 1-alpha
   # we derive alpha = 1 - prob
   alpha = 1 - prob
   prb = 1.0 - alpha/2
   tval = t.ppf(prb, fitobj.dof)

   C = fitobj.covar
   n = len(fitobj.params)              # Number of parameters from covariance matrix
   p = fitobj.params
   N = len(x)
   if abswei:
      covscale = 1.0  # Do not apply correction with red. chi^2
   else:
      covscale = fitobj.rchi2_min
   df2 = numpy.zeros(N)
   for j in range(n):
      for k in range(n):
         df2 += dfdp[j]*dfdp[k]*C[j,k]
   if prediction:
      df = numpy.sqrt(err*err+covscale*df2)
   else:
      df = numpy.sqrt(covscale*df2)
   y = f(p, x)
   delta = tval * df
   upperband = y + delta
   lowerband = y - delta
   return y, upperband, lowerband

def confidence_band(x, dfdp, confprob, fitobj, f, err=None, abswei=False):
   #----------------------------------------------------------
   # Given a value for x, calculate the error df in y = model(p,x)
   # This function returns for each x in a NumPy array, the
   # upper and lower value of the confidence interval.
   # The arrays with limits are returned and can be used to
   # plot confidence bands.
   #
   #
   # Input:
   #
   # x        NumPy array with values for which you want
   #          the confidence interval.
   #
   # dfdp     A list with derivatives. There are as many entries in
   #          this list as there are parameters in your model.
   #
   # confprob Confidence probability in percent (e.g. 90% or 95%).
   #          From this number we derive the confidence level
   #          (e.g. 0.05). The Confidence Band
   #          is a 100*(1-alpha)% band. This implies
   #          that for a given value of x the probability that
   #          the 'true' value of f(p,x) falls within these limits is
   #          100*(1-alpha)%.
   #
   # fitobj   The Fitter object from a fit with kmpfit
   #
   # f        A function that returns a value y = f(p,x)
   #          p are the best-fit parameters and x is a NumPy array
   #          with values of x for which you want the confidence interval.
   #
   # abswei   Are the weights absolute? For absolute weights we take
   #          unscaled covariance matrix elements in our calculations.
   #          For unit weighting (i.e. unweighted) and relative
   #          weighting, we scale the covariance matrix elements with
   #          the value of the reduced chi squared.
   #
   # Returns:
   #
   # y          The model values at x: y = f(p,x)
   # upperband  The upper confidence limits
   # lowerband  The lower confidence limits
   #
   # Note:
   #
   # If parameters were fixed in the fit, the corresponding
   # error is 0 and there is no contribution to the confidence
   # interval.
   #----------------------------------------------------------
   return confpred_band(x, dfdp, confprob, fitobj, f, prediction=False, err=err, abswei=abswei)

def prediction_band(x, dfdp, predprob, fitobj, f, err=None, abswei=False):
   #----------------------------------------------------------
   # Given a value for x, calculate the error df in y = model(p,x)
   # This function returns for each x in a NumPy array, the
   # upper and lower value of the prediction interval.
   # The arrays with limits are returned and can be used to
   # plot confidence bands.
   #
   #
   # Input:
   #
   # x        NumPy array with values for which you want
   #          the prediction interval.
   #
   # dfdp     A list with derivatives. There are as many entries in
   #          this list as there are parameters in your model.
   #
   # predprob Prediction probability in percent (e.g. 0.9 or 0.95).
   #          From this number we derive the prediction level
   #          (e.g. 0.05). The Prediction Band
   #          is a 100*(1-alpha)% band. This implies
   #          that values of one or more future observations from
   #          the same population from which a given data set was sampled,
   #          will fall in this band with a probability of 100*(1-alpha)%
   #
   # fitobj   The Fitter object from a fit with kmpfit
   #
   # f        A function that returns a value y = f(p,x)
   #          p are the best-fit parameters and x is a NumPy array
   #          with values of x for which you want the confidence interval.
   #
   # abswei   Are the weights absolute? For absolute weights we take
   #          unscaled covariance matrix elements in our calculations.
   #          For unit weighting (i.e. unweighted) and relative
   #          weighting, we scale the covariance matrix elements with
   #          the value of the reduced chi squared.
   #
   # Returns:
   #
   # y          The model values at x: y = f(p,x)
   # upperband  The upper prediction limits
   # lowerband  The lower prediction limits
   #
   # Note:
   #
   # If parameters were fixed in the fit, the corresponding
   # error is 0 and there is no contribution to the prediction
   # interval.
   #----------------------------------------------------------
   return confpred_band(x, dfdp, predprob, fitobj, f,
                        prediction=True, err=err, abswei=abswei)

def my_model(p, x):
   #-----------------------------------------------------------------------
   # This describes the model and its parameters for which we want to find
   # the best fit. 'p' is a sequence of parameters (array/list/tuple).
   #-----------------------------------------------------------------------
   A, mu, sigma, zerolev = p
   return( A * numpy.exp(-(x-mu)*(x-mu)/(2.0*sigma*sigma)) + zerolev )

def my_residuals(p, data):
   #-----------------------------------------------------------------------
   # This function is the function called by the fit routine in kmpfit
   # It returns a weighted residual. De fit routine calculates the
   # square of these values.
   #-----------------------------------------------------------------------
   x, y, err = data
   return (y-my_model(p,x)) / err

def my_derivs(p, data, dflags):
   #-----------------------------------------------------------------------
   # This function is used by the fit routine to find the values for
   # the explicit partial derivatives. Argument 'dflags' is a list
   # with booleans. If an element is True then an explicit partial
   # derivative is required.
   #-----------------------------------------------------------------------
   x, y, err = data    # y is dummy here
   A, mu, sigma, zerolev = p
   pderiv = numpy.zeros([len(p), len(x)])  # You need to create the required array
   sig2 = sigma * sigma
   sig3 = sig2 * sigma
   xmu  = x-mu
   xmu2 = xmu**2
   expo = numpy.exp(-xmu2/(2.0*sig2))
   fx = A * expo
   for i, flag in enumerate(dflags):
      if flag:
         if i == 0:
            pderiv[0] = expo
         elif i == 1:
            pderiv[1] = fx * xmu/(sig2)
         elif i == 2:
            pderiv[2] = fx * xmu2/(sig3)
         elif i == 3:
            pderiv[3] = 1.0
   return pderiv/-err

# Artificial data
N = 50
x = numpy.linspace(-5, 10, N)
truepars = [10.0, 5.0, 1.0, 0.0]
p0 = [9, 4.5, 0.8, 0]
rms_data = 0.8
rms_err = 0.1
y = my_model(truepars, x) + numpy.random.normal(0.0, rms_data, N)
err = numpy.random.normal(0.6, rms_err, N)
#err = err*0 + 1

# The fit
fitobj = kmpfit.Fitter(residuals=my_residuals, deriv=my_derivs, data=(x, y, err))
try:
   fitobj.fit(params0=p0)
except Exception as mes:
   print("Something wrong with fit: ", mes)
   raise SystemExit

print("\n\n======== Results kmpfit with explicit partial derivatives =========")
print("Params:        ", fitobj.params)
print("Errors from covariance matrix         : ", fitobj.xerror)
print("Uncertainties assuming reduced Chi^2=1: ", fitobj.stderr)
print("Chi^2 min:     ", fitobj.chi2_min)
print("Reduced Chi^2: ", fitobj.rchi2_min)
print("Iterations:    ", fitobj.niter)
print("Function ev:   ", fitobj.nfev)
print("Status:        ", fitobj.status)
print("Status Message:", fitobj.message)
print("Covariance:\n", fitobj.covar)

# Re-use my_derivs() but rescale derivatives back again with -err
dervs = my_derivs(fitobj.params, (x,y,err), (True,True,True,True))*-err

dfdp = [dervs[0], dervs[1], dervs[2], dervs[3]]
confprob = 0.95
ydummy, upperband, lowerband = confidence_band(x, dfdp, confprob, fitobj, my_model)
verts_conf = list(zip(x, lowerband)) + list(zip(x[::-1], upperband[::-1]))

predprob = 0.90
ydummy, upperband, lowerband = prediction_band(x, dfdp, predprob, fitobj, my_model,
                               err=err, abswei=False)
verts_pred = list(zip(x, lowerband)) + list(zip(x[::-1], upperband[::-1]))

# Plot the result
rc('font', size=9)
rc('legend', fontsize=8)
fig = figure()
frame = fig.add_subplot(1,1,1)
X = numpy.linspace(x.min(), x.max(), 100)
frame.errorbar(x, y, yerr=err, fmt='go', alpha=0.7, label="Noisy data")
frame.plot(X, my_model(truepars,X), 'r', label="True data")
frame.plot(X, my_model(fitobj.params,X), 'b', lw=2, label="Fit with kmpfit")
poly = Polygon(verts_conf, closed=True, fc='g', ec='g', alpha=0.3,
               label="CI (%g)"%confprob)
frame.add_patch(poly)
poly = Polygon(verts_pred, closed=True, fc='r', ec='r', alpha=0.3,
               label="PI (%g)"%predprob)
frame.add_patch(poly)
frame.set_xlabel("X")
frame.set_ylabel("Measurement data")
frame.set_title("Confidence- and prediction bands for Gaussian model",
                fontsize=10)
delta = (x.max()-x.min())/10.0
frame.set_xlim(x.min()-delta, x.max()+delta)
frame.grid(True)

# Check prediction intervals
"""
for i in range(500):
   y = my_model(truepars, x) + numpy.random.normal(0.0, rms_data, N)
   err = numpy.random.normal(0.0, rms_err, N)
   #frame.plot(x,y,'o')
   frame.errorbar(x, y, yerr=err, fmt='o')
"""
# A nice background for the entire plot
from matplotlib.cm import copper
frame.imshow([[0, 0],[1,1]], interpolation='bicubic', cmap=copper,
             vmin=-0.5, vmax=0.5,
             extent=(frame.get_xlim()[0], frame.get_xlim()[1],
                     frame.get_ylim()[0], frame.get_ylim()[1]),
             alpha=1)

leg = frame.legend(loc=2)
show()

max_stars_repo_name: kapteyn-astro/kapteyn

max_issues_repo_name: kapteyn-astro/kapteyn

max_forks_repo_name: kapteyn-astro/kapteyn

reason: import numpy,from scipy


text:
import itertools
import subprocess
from io import StringIO
import shutil
import mdtraj as md
from simtk import openmm as mm
import mdtraj.utils as mdtrajutils
import numpy as np
import openmmtools
import pandas as pd
from openmmtools.mcmc import WeightedMove, LangevinSplittingDynamicsMove, SequenceMove, \
    MCDisplacementMove, MCRotationMove, GHMCMove
from simtk import unit
from simtk.openmm import app

from rlmm.environment.systemloader import working_directory
from rlmm.utils.config import Config

def get_mdtraj_box(boxvec=None, a=None, b=None, c=None, iterset=-1):
    if boxvec is not None:
        a, b, c = boxvec
    elif None in [a, b, c]:
        exit()

    a, b, c = a.value_in_unit(unit.angstrom), b.value_in_unit(unit.angstrom), c.value_in_unit(unit.angstrom)
    a, b, c = np.array(a), np.array(b), np.array(c)
    a, b, c, alpha, beta, gamma = md.utils.unitcell.box_vectors_to_lengths_and_angles(a, b, c)
    if iterset is -1:
        return a, b, c, alpha, beta, gamma
    else:
        return [[a, b, c]] * iterset, [[alpha, beta, gamma]] * iterset

def get_coordinates_samplers(topology: app.topology.Topology, sampler_state: openmmtools.states.SamplerState, explicit: bool):
    return get_coordinates(topology, sampler_state.positions, sampler_state.box_vectors, sampler_state.n_particles,
                           explicit)

def get_coordinates(topology: app.topology.Topology,
                    coords,
                    boxvecs,
                    n_particles: int,
                    explicit: bool):
    coords = np.array(coords.value_in_unit(unit.angstrom))

    if explicit:
        trajectory_positions = coords.reshape(
            (1, n_particles, 3))
        lengths, angles = get_mdtraj_box(boxvec=boxvecs, iterset=1)

        traj = md.Trajectory(trajectory_positions, md.Topology.from_openmm(topology),
                             unitcell_lengths=lengths, unitcell_angles=angles)

        traj = traj.image_molecules(inplace=False)
        coords = traj.xyz.reshape((traj.n_atoms, 3))

    return coords

def get_pdb(topology, coords, file_name=None):
    if file_name is None:
        output = StringIO()
    else:
        output = open(file_name, 'w')
    app.PDBFile.writeFile(topology,
                          coords,
                          file=output, keepIds=False)
    if file_name is None:
        return output.getvalue()
    else:
        output.close()
        return True

def get_ligand_ids(topology):
    return md.Topology.from_openmm(topology).select("resn UNL")

def get_protein_ids(topology):
    return md.Topology.from_openmm(topology).select("protein")

def get_backbone_ids(topology):
    return md.Topology.from_openmm(topology).select("protein and backbone")

def get_selection_ids(topology, sele):
    return md.Topology.from_openmm(topology).select(sele)

def get_protein_restraint_force(topology, positions, explicit, K=5.0):
    if explicit:
        energy_expression = '(k_restr/2)*periodicdistance(x, y, z, x0, y0, z0)^2' # periodic distance
    else:
        energy_expression = '(k_restr/2)*((x-x0)^2 + (y-y0)^2 + (z-z0)^2)' # non-periodic distance

    force = mm.CustomExternalForce(energy_expression)
    force.addGlobalParameter("k_restr", K)
    force.addPerParticleParameter("x0")
    force.addPerParticleParameter("y0")
    force.addPerParticleParameter("z0")
    positions_ = positions.value_in_unit(unit.nanometer)
    for i, atom_id in enumerate(md.Topology.from_openmm(topology).select("protein and (not type H)")):
        pos = positions_[atom_id]
        pops = mm.Vec3(pos[0], pos[1], pos[2])
        _ = force.addParticle(int(atom_id), pops)
    return force

def get_backbone_restraint_force(topology, positions, explicit, K=5.0):
    if explicit:
        energy_expression = '(k_restr/2)*periodicdistance(x, y, z, x0, y0, z0)^2' # periodic distance
    else:
        energy_expression = '(k_restr/2)*((x-x0)^2 + (y-y0)^2 + (z-z0)^2)' # non-periodic distance

    force = mm.CustomExternalForce(energy_expression)
    force.addGlobalParameter("k_restr", K)
    force.addPerParticleParameter("x0")
    force.addPerParticleParameter("y0")
    force.addPerParticleParameter("z0")
    positions_ = positions.value_in_unit(unit.nanometer)
    for i, atom_id in enumerate(get_backbone_ids(topology)):
        pos = positions_[atom_id]
        pops = mm.Vec3(pos[0], pos[1], pos[2])
        _ = force.addParticle(int(atom_id), pops)
    return force

def get_ligand_restraint_force(topology, positions, explicit, K=5.0):
    if explicit:
        energy_expression = '(k_restr/2)*periodicdistance(x, y, z, x0, y0, z0)^2' # periodic distance
    else:
        energy_expression = '(k_restr/2)*((x-x0)^2 + (y-y0)^2 + (z-z0)^2)' # non-periodic distance

    force = mm.CustomExternalForce(energy_expression)
    force.addGlobalParameter("k_restr", K)
    force.addPerParticleParameter("x0")
    force.addPerParticleParameter("y0")
    force.addPerParticleParameter("z0")
    positions_ = positions.value_in_unit(unit.nanometer)
    for i, atom_id in enumerate(md.Topology.from_openmm(topology).select("resn UNL and (not type H)")):
        pos = positions_[atom_id]
        pops = mm.Vec3(pos[0], pos[1], pos[2])
        _ = force.addParticle(int(atom_id), pops)
    return force

def get_pocket_residues(traj):
    traj = traj.atom_slice(traj.topology.select("protein or resn UNL"))
    resn = len(list(traj.topology.residues))
    group_1 = list(range(resn-1))
    group_2 = [resn - 1]
    pairs = list(itertools.product(group_1, group_2))
    res, pairs = md.compute_contacts(traj, pairs)
    pocket_resids = list(np.where(res[0] <= 0.5)[0] + 1)
    pocket_resids = ["resid {}".format(id) for id in pocket_resids]
    pocket_resids = " or ".join(pocket_resids)
    return pocket_resids

def detect_ligand_flyaway(traj, pocket_resids, eps=0.2, return_difference=False):
    group_1 = list(traj.topology.select(pocket_resids))
    group_2 = list(traj.topology.select("resn UNL"))
    pairs = list(itertools.product(group_1, group_2))
    res = md.compute_distances(traj, pairs)
    distances = np.quantile(res, 0.99, axis=1)
    difference = np.abs(np.mean(distances[:int(distances.shape[0] * 0.05)]) - np.mean(distances[int(distances.shape[0]* 0.95):]))
    return difference >= eps if not return_difference else (difference >= eps, difference)

def run_amber_mmgbsa(logger, explicit, tempdir, run_decomp=False):

    with logger('run_amber_mmgbsa') as logger:
        logger.log("PASSING")
        return
        complex_prmtop = f"com.prmtop"
        traj = "traj.dcd"
        with working_directory(tempdir):
            if explicit:
                with open("cpptraj_input.txt", 'w') as f:
                    f.write("strip :WAT parmout stripped.prmtop outprefix traj.dcd nobox\n" +
                            "trajout test2.dcd\n" +
                            "run\n")
                subprocess.run(['cpptraj', '-p', complex_prmtop, '-y', traj, '-i', 'cpptraj_input.txt'],
                                      check=True, capture_output=True)
                complex_prmtop = "stripped.prmtop"
                traj = "test2.dcd"
            else:
                shutil.copyfile(com)

            try:
                subprocess.run(['ante-MMPBSA.py',
                                       '-p', complex_prmtop,
                                       '-l', 'noslig.prmtop',
                                       '-r', 'nosapo.prmtop',
                                       '-n', ':UNL'], check=True, capture_output=True)
            except subprocess.CalledProcessError as e:
                logger.error(e, e.stderr.decode("UTF-8"), e.stdout.decode("UTF-8"))

            with open("mmpbsa_input.txt", 'w') as f:
                f.write(
                    '&general\ninterval=5,\nverbose=3, keep_files=0, strip_mask=":WAT:CL:CIO:CS:IB:K:LI:MG:NA:RB:HOH",\n/\n&gb\nigb=5, saltcon=0.1000,\n/\n'
                )
                if run_decomp:
                    f.write('&decomp\nidecomp=1,csv_format=1\n/\n')

            logger.log("Running amber MMPBSA.py, might take awhile...")
            # process = subprocess.Popen(['MMPBSA.py', '-y', traj,
            #                        '-i', 'mmpbsa_input.txt',
            #                        '-cp', complex_prmtop,
            #                        '-rp', 'nosapo.prmtop',
            #                        '-lp', 'noslig.prmtop'], stdout=subprocess.DEVNULL, stderr=subprocess.DEVNULL)
            #
            subprocess.run(['MMPBSA.py', '-y', traj,
                                   '-i', 'mmpbsa_input.txt',
                                   '-cp', complex_prmtop,
                                   '-rp', 'nosapo.prmtop',
                                   '-lp', 'noslig.prmtop'], capture_output=True, check=True)
            #
            # decomp = self.decomp_to_csv('FINAL_DECOMP_MMPBSA.dat', 'decomp.csv')
            results = results_to_csv('FINAL_RESULTS_MMPBSA.dat', 'result.csv')
            logger.log(results.iloc[-1])

def decomp_to_csv(decomp_filename, csv_filename):
    with open(decomp_filename, 'r') as f:
        header = ''
        while 'Resid 1' not in header:
            header = f.readline()
        h1 = header.strip().split(',')
        h2 = f.readline()
        h2 = h2.strip().split(',')
        resid = h1[:2]
        h1, h2 = list(dict.fromkeys(h1[2:])), list(dict.fromkeys(h2))
        h1.remove('')
        h2.remove('')
        header = ','.join(resid + list(map(' '.join, list(itertools.product(h1, h2))))) + '\n'
        with open(csv_filename, 'w') as csvfile:
            csvfile.write(header)
            csvfile.writelines(f.readlines())
    return pd.read_csv(csv_filename)

def results_to_csv(results_filename, csv_filename):
    with open(results_filename, 'r') as f:
        header = ''
        while 'Differences (Complex - Receptor - Ligand)' not in header:
            header = f.readline()
        header = f.readline()
        header = [h.strip() for h in header.strip().split('  ') if h]
        header = ','.join(header) + '\n'
        with open(csv_filename, 'w') as csvfile:
            csvfile.write(header)
            f.readline()
            while True:
                l = f.readline()
                if '---' not in l:
                    l = [c.strip() for c in l.split('   ') if c]
                    l = ','.join(l) + '\n'
                    csvfile.write(l)
                if 'TOTAL' in l:
                    break
    return pd.read_csv(csv_filename)

def prepare_mcmc(topology, config):
    sequence_move = LangevinSplittingDynamicsMove(
        timestep=config.parameters.integrator_params['timestep'],
        n_steps=config.n_steps,
        collision_rate=config.parameters.integrator_params['collision_rate'],
        reassign_velocities=True,
        n_restart_attempts=6,
        constraint_tolerance=config.parameters.integrator_setConstraintTolerance)

    if config.hybrid:
        atoms = md.Topology.from_openmm(topology).select("resn UNL")
        rot_move = MCRotationMove(atom_subset=atoms)
        trans_move = MCDisplacementMove(atom_subset=atoms,displacement_sigma=config.displacement_sigma * unit.angstrom)
        sequence_move = SequenceMove([trans_move, rot_move, sequence_move])

    return sequence_move

class SystemParams(Config):
    def __init__(self, config_dict):
        # noinspection PyUnresolvedReferences
        import simtk.openmm as mm
        self.createSystem = None
        self.platform = None
        self.minMaxIters = None
        self.integrator_setConstraintTolerance = None
        self.platform_config = None
        self.integrator_params = None
        for k, v in config_dict.items():
            if k != "platform_config" and isinstance(v, dict):
                for k_, v_ in v.items():
                    try:
                        exec('v[k_] = ' + v_)
                    except TypeError:
                        v[k_] = v_
            else:
                exec('config_dict[k] = ' + str(v))
        self.__dict__.update(config_dict)

ext: py

max_stars_repo_name: dallonasnes/RLMM

max_issues_repo_name: dallonasnes/RLMM

max_forks_repo_name: dallonasnes/RLMM

reason: import numpy


text:
import cv2
import numpy as np
from matplotlib import pyplot as plt

img = cv2.imread('img/test.png',0)
f = np.fft.fft2(img)
fshift = np.fft.fftshift(f)
magnitude_spectrum = 20*np.log(np.abs(fshift))

plt.subplot(121),plt.imshow(img, cmap = 'gray')
plt.title('Input Image'), plt.xticks([]), plt.yticks([])
plt.subplot(122),plt.imshow(magnitude_spectrum, cmap = 'gray')
plt.title('Magnitude Spectrum'), plt.xticks([]), plt.yticks([])
plt.show()

max_stars_repo_name: ThomasDerZweifler/pyPro

max_issues_repo_name: ThomasDerZweifler/pyPro

max_forks_repo_name: ThomasDerZweifler/pyPro


text:
from itertools import groupby
from queue import Queue
import time
import sys

import pandas as pd
from enum import Enum
from collections import OrderedDict

import utils
from prog import re_allocate
import psutil
import numpy as np
import managerOperations
from msgOperations import MSGOperations
import sys_admin
from threading import Thread
import trace_generator

MAXIMUM_PARALLEL = 2 #

MONITOR_GAP = 20
NUM_OF_GPUs_PER_NODE = 8 #

class JobNodeStatus(Enum):
    JOBFAILED=0
    JOBOUT=1
    NODEIN=2
    NODEOUT=3

class Manager:
    def __init__(self, max_parallel = MAXIMUM_PARALLEL, monitor_gap = MONITOR_GAP):
        self.max_parallel = max_parallel
        self.monitor_gap = monitor_gap
        self.create_working_directory()
        self.buffer = Queue()

        # create server ready to recv data
        self.run_msg_server()

    current_map = pd.DataFrame()
    # A dictionary for global job information
    job_info_dict = {}

    # create working directory
    def create_working_directory(self):
        managerOperations.create_working_directory()

    # Network Operations
    def create_msg_client(self, address, port):
        """Create a message client by client for reporting training throughput/speed.

        Args:
            func (func pointer): Client need to offer a function to get the real time training speed
        """
        return MSGOperations().create_msg_client(address, port)

    def create_msg_server(self): # pass dynamic update data function into
        """Create a message server for receive throughput report by client

        Args:
            func (function pointer): a function for deciding when to trigger update job data
        """
        managerOperations.create_msg_server()

    # life cycle functions
    def manager_start(self):
        sys_nodes=sys_admin.get_cluster_nodes()
        if len(sys_nodes) == 0 or (sys_nodes is None):
            return

        if utils.get_job_queue_len() == 0:
            return

        # fetch job from DB
        starting_jobs_num = min(MAXIMUM_PARALLEL, utils.get_job_queue_len())
        for i in range(starting_jobs_num):
            job_string = utils.get_a_job_from_DB()
            jobdetail = utils.parser_job_string_2_job_item(job_string)
            self.job_info_dict[jobdetail.GUID] = jobdetail

        # build inital map
        jobnames = self.job_info_dict.keys()
        data = np.zeros((len(jobnames), len(sys_nodes)), dtype=int) # initial fake data
        initialMap = pd.DataFrame(data=data, index=jobnames, columns=sys_nodes)

        # Feed the optimizer and get the optimized newmap
        mins, maxs, Ns, Os, res_ups, res_dws = utils.get_optimizer_parameters_by_job_dict(self.job_info_dict)

        tmpGRB, new_data, tmpRate, tmpCost = re_allocate(cmap=initialMap, jmin=mins, jmax=maxs,
                                                            Ns=Ns,Os=Os, Tfwd=10, res_up=res_ups, res_dw = res_dws,
                                                            time_limit=10)
        print(new_data)

        new_map = pd.DataFrame(data=new_data, index=jobnames, columns=sys_nodes)

        managerOperations.adjust_nodes_by_map(new_map=new_map, old_map=initialMap, job_info_dict = self.job_info_dict)
        self.current_map = new_map

        print("=======================")
        print("===manager start end===")
        print("=======================")

    def scheduler_job_change(self, GUIDs):

        # 1. Detect job leave
        # drop the job from map and job_dict
        for GUID in GUIDs:
            print("GUID: ", GUID)
            print("GUID Type: ", type(GUID))
            # 1.1 Drop job from current map
            self.current_map.drop([GUID], inplace=True)
            # 1.2 Delete jobInfo from jobInfo dictionary
            self.job_info_dict.pop(GUID)

        # 2. fetch jobs to max parallel(batch manner)
        print("job leaving after and fetch before job_info_dict: ", self.job_info_dict)
        print("fetching job to max parallel")

        job_string_list = []
        lacking_len = self.max_parallel - len(self.job_info_dict)
        print(self.max_parallel)
        print("lacking_len: ", lacking_len)

        left_jobs_in_db = utils.get_job_queue_len()
        print("There are %d jobs left in database" % left_jobs_in_db)

        job_num = min(lacking_len, left_jobs_in_db)
        print("fetched valid job number:", job_num)

        for i in range(job_num):
            jobstring = utils.get_a_job_from_DB()
            print("fetch new job: ", jobstring)
            job_string_list.append(jobstring)
            jobdetail = utils.parser_job_string_2_job_item(jobstring)
            newrow = pd.DataFrame([np.zeros(len(self.current_map.columns), dtype=int)], index=[jobdetail.GUID] ,columns=self.current_map.columns)

            # add new job to map
            self.current_map = self.current_map.append(newrow)
            # add new job to dict
            self.job_info_dict[jobdetail.GUID] = jobdetail

        print("After fetching jobs get the map for adjustment")
        print(self.current_map)

        # update buffer info to job_info_dict
        #
        self.update_job_data_on_events(self.buffer, event_type="job event")

        # get parameters with latest `job_info_dict`
        mins, maxs, Ns, Os, res_ups, res_dws = utils.get_optimizer_parameters_by_job_dict(self.job_info_dict)

        tmpGRB, new_data, tmpRate, tmpCost = re_allocate(cmap=self.current_map, jmin=mins, jmax=maxs,
                                Ns=Ns,Os=Os, Tfwd=10, res_up=res_ups, res_dw = res_dws, time_limit=10)

        new_map = pd.DataFrame(data=new_data, index=self.current_map.index, columns=self.current_map.columns)

        managerOperations.adjust_nodes_by_map(new_map, self.current_map, self.job_info_dict)

        # update current_map
        self.current_map = new_map

        # print current map after allocation
        print("After re-allocation jobs get the new map")
        print(self.current_map)

    def scheduler_nodes_change(self, flag, nodes):
        print("node change called")

        # Event driven update data here
        # Use buffer data and get ns os and res_up and res_dw
        #
        self.update_job_data_on_events(self.buffer, event_type= "node event")

        # validate nodes name before operations
        if sys_admin.is_nodes_belong_to_avaliable_nodes(nodes) == False:
            print("error: nodes out of avaliable nodes range")
            return

        # drop or add columns for nodes in cmap
        old_map = self.current_map
        mins, maxs, Ns, Os, res_ups, res_dws = utils.get_optimizer_parameters_by_job_dict(self.job_info_dict)

        print("Ns Os res_ups res_dw")
        print(Ns, Os, res_ups, res_dws)

        if flag == JobNodeStatus.NODEIN:
            print("node in :", nodes)
            print("old map", self.current_map)
            for node in nodes:
               self.current_map.insert(self.current_map.shape[1], node, 0) # dataframe add one new column
            tmpGRB, new_data, tmpRate, tmpCost = re_allocate(cmap=self.current_map, jmin=mins, jmax=maxs,
                                                                Ns=Ns,Os=Os, Tfwd=10, res_up=res_ups, res_dw = res_dws, time_limit=10)
            new_map = pd.DataFrame(data=new_data, index=self.current_map.index, columns=self.current_map.columns)
        else:
            print("node leave ", nodes)
            print("old map", self.current_map)
            for node in nodes:
                tmp_map = self.current_map.drop(labels=node, axis=1) # dataframe delete one column

            tmpGRB, new_data, tmpRate, tmpCost = re_allocate(cmap=tmp_map, jmin=mins, jmax=maxs,
                                                Ns=Ns,Os=Os, Tfwd=10, res_up=res_ups, res_dw = res_dws, time_limit=10)
            new_map = pd.DataFrame(data=new_data, index=tmp_map.index, columns=tmp_map.columns)
            print("new map", new_map)

        managerOperations.adjust_nodes_by_map(new_map, old_map, self.job_info_dict)

        # update current_map
        self.current_map = new_map

    # for future usage
    def _terminate_manager(self):
        # kill all horovodrun process in dictionary
        for job in self.job_info_dict:
            p = psutil.Process(self.job_info_dict[job].pid)
            p.terminate() # kill all the running process for this job

        # clean all host file and discovery file
        for job in self.job_info_dict:
            managerOperations.del_host_files(jobname=job)

    def monitor_hvd_processes(self):
        while True:
            now = time.time()
            print("============= %f : monitor hvd process report * Head ===============" % now)

            time.sleep(self.monitor_gap)
            jobnames = []
            print(self.current_map)
            for jobname in self.job_info_dict.keys():
                pid = self.job_info_dict[jobname].pid
                print("dict stored hvd job - pid: %s - %s" % (jobname, pid))

                if not psutil.pid_exists(pid):

                    print("Hey! we found hvdrun job: %s - process: %s not existed anymore, so we start fetch new jobs" % (jobname, pid))
                    jobnames.append(jobname)

            if jobnames:
                print("jobnames", jobnames)
                self.scheduler_job_change(GUIDs=jobnames)

            print("============= %f : monitor hvd process report * foot ===============" % now)

    def merge_ordered_NO_2_itemNO(self, job_N, job_O, N, O):
        # Switch job_N and job_O to a dict
        item_no_dict = {}
        for i in range(len(job_N)):
            item_no_dict[job_N[i]] = job_O[i]

        # insert N and O to the tmp dict
        # Key same will replace, otherwise will insert
        for i in range(len(N)):
            print("insert key Ni", N[i])
            print("insert val Oi", O[i])
            item_no_dict[N[i]] = O[i]

        # Sort the merged dict
        ordered_job_dict = OrderedDict(sorted(item_no_dict.items()))

        # convert back to N and O (ordered)
        Ns = list(ordered_job_dict.keys())
        Os = list(ordered_job_dict.values())

        # print res_N and res_O
        print("Ns: ", Ns)
        print("Os: ", Os)

        return Ns, Os

    def update_scaling_and_cost_data_2_jobInfoDict(self, jobname, N=None, O=None, res_up = None, res_down = None):
        """
        This function only for dynamic update job_info_dict job N O and resup and down information
        """

        print("========update_scaling_and_cost_data_2_jobInfoDict called==========")
        print("N is: ", N)
        print("O is: ", O)
        print("res_up is", res_up)
        print("res_dw is", res_down)

        if self.job_info_dict == None or len(self.job_info_dict) == 0:
            return

        if jobname == None or len(jobname) == 0:
            return

        job_item = self.job_info_dict[jobname] # assume we can find this job info by name

        if N != None or O != None:
            if job_item.N == None or job_item.O == None:
                job_item.N = N
                job_item.O = O
            else:
                # if not None then get the merged N and O
                res_N, res_O = self.merge_ordered_NO_2_itemNO(job_N=job_item.N, job_O=job_item.O, N=N, O=O)

                # update job item
                job_item.N = res_N
                job_item.O = res_O

        if res_up != None:
            job_item.res_up = res_up

        if res_down != None:
            job_item.res_down = res_down

    def _get_group_dict(self, buffer):
        # get job msg_item dict
        group_dict = {}
        for key in buffer:
            hostname = utils.get_host_name_by_address(key)
            jobname = utils.get_jobname_by_hostname(hostname, self.current_map)

            if jobname == "":
                continue

            msg_list = buffer[key]
            msg_items = []
            for msg in msg_list:
                msg_item = utils.parser_udp_message(msg)
                msg_items.append(msg_item)
            group_dict[jobname] = msg_items

        return group_dict

    def _cal_res_up_res_dw(self, job_items):
        res_up, res_dw = None, None
        if len(job_items) > 2:
            print("update res_up and res_dw")

            for i in range(len(job_items)-1,-1,-1): # reverse order find rank difference

                if job_items[i].rank_size > job_items[i-1].rank_size and job_items[i-1].rank_size == job_items[i-2].rank_size:
                    print("================ get the res_up cost ===================")
                    res_up = (job_items[i].time - job_items[i-1].time) - (job_items[i-1].time - job_items[i-2].time)
                    print("==resup==", res_up)

                elif job_items[i].rank_size < job_items[i-1].rank_size and job_items[i-1].rank_size == job_items[i-2].rank_size:
                    print("================ get the res_down cost =================")
                    print(type(job_items[i].time))
                    print(job_items[i].time)
                    res_dw = (job_items[i].time - job_items[i-1].time) - (job_items[i-1].time - job_items[i-2].time)
                    print("==resdw==", res_dw)
        return res_up, res_dw

    def update_job_data_on_events(self, buffer, event_type):
        ''' job change or node change trigger updating data for re-allocation'''
        if len(buffer) == 0:
            print("No valid information and skip update process")
            return

        group_dict = self._get_group_dict(buffer)

        if len(group_dict) == 0:
            print("Jobs in current_map has no any training msg in buffer, probably means all jobs are new fetched(not started yet)," \
            " do not need to use historial information to do the update, so return the update function here")
            return

        # for each job get the cost and thrpt info
        for jobname in group_dict:
            job_items = group_dict[jobname]

            # we will assume the id is right for now
            # job_items.sort(key=lambda x: x.id)

            res_up, res_dw = self._cal_res_up_res_dw(job_items)

            # TODO: the cal throughput part need to change(use latest throughput)
            N = []
            O = []
            for key, group in groupby(job_items, lambda x: x.rank_size): # key: different rank size for job - group: items of this job with this ranksize
                node_num = int(key)/NUM_OF_GPUs_PER_NODE
                N.append(int(node_num))
                group_list = list(group)

                thrputs = []
                for i in range(0, len(group_list) - 1):
                    msg_time_gap = float(group_list[i + 1].time) - float(group_list[i].time)
                    thrput = float(group_list[i].credit) / msg_time_gap
                    thrputs.append(thrput)
                avg_thrput = thrputs[-1] # use the last thrput as the current thrput

                O.append(avg_thrput)

        # Update collect info to JobInfoDict
        self.update_scaling_and_cost_data_2_jobInfoDict(jobname=jobname, N=N, O=O, res_up=res_up, res_down=res_dw)

    def run_msg_server(self):
        # run server daemon
        print("run server and update data")
        mserver=MSGOperations()
        self.buffer = mserver.buffer

        p_server = Thread(target=mserver.create_msg_server) # keep updating buffer data
        p_server.start()

    def is_first_round(self, counters):
        flag = True
        for coun in list(counters.values()):
            if coun != 0:
                flag = False
        return flag

    # node come and leave
    def events_launcher(self):
        self.manager_start()

        # create events source
        cluster_nodes = sys_admin.get_cluster_nodes()
        trace = trace_generator.synthetic_trace(nodes=cluster_nodes, nf=20000)

        # write events trace to log
        with open("events.log", 'w') as f:
            f.write(str(trace))

        # record the status of nodes in arrary for controlling events
        counters = {}
        flags = {}
        for node in cluster_nodes:
            counters[node] = 0
            flags[node] = False

        start_time = time.time()

        # events tick driver
        while(True):
            time.sleep(0.1) # check each 0.1 second(could change)

            passing_time = time.time() - start_time

            coming_nodes = [] # coming nodes in one time check
            leaving_nodes = [] # leaving nodes in one time check

            for node in cluster_nodes:
                timestamps_tuple = trace[node] # list of tuple for node
                current_time_tuple = timestamps_tuple[counters[node]]

                if passing_time > current_time_tuple[0] and flags[node] == False:
                    print("node: " + node + " in") # trigger node in event
                    coming_nodes.append(node)
                    flags[node] = True

                if passing_time > current_time_tuple[1] and flags[node] == True:
                    print("node: " + node + " leave") # trigger node leave event
                    leaving_nodes.append(node)
                    flags[node] = False
                    counters[node] = counters[node] + 1

            # when node counter == 0 that is the start, just skip the start phase
            if coming_nodes and self.is_first_round(counters) == False:
                self.scheduler_nodes_change(JobNodeStatus.NODEIN, coming_nodes)

            if leaving_nodes:
                self.scheduler_nodes_change(JobNodeStatus.NODEOUT, leaving_nodes)

def main():

    # create manager
    # basic settings and msg server
    m = Manager()

    # 2. start jobs and run as events come
    p_events = Thread(target=m.events_launcher)
    p_events.start()

    # 3. process monitor job pid
    p_monitor = Thread(target=m.monitor_hvd_processes)
    p_monitor.start()

    # for local testing purpose

if __name__ == "__main__":
    main()

max_stars_repo_name: BFTrainer/BFTrainer

max_issues_repo_name: BFTrainer/BFTrainer

max_forks_repo_name: BFTrainer/BFTrainer


text:
"""predictions of users based on tweet embeddings"""
import numpy as np
from sklearn.linear_model import LogisticRegression
from .twitter import vectorize_tweet
from .models import User

def predict_user(user0_name, user1_name, hypo_tweet_text):
    """
    Determine and return which user is more likely to say a hypothetical tweet.

    example run: predicts_user("elonmusk", "jackblakc",
                               "school of rock really rocks")
    returns 0 (user0_name=elonmusk) or 1 (user1name=jackblack)
    """

    # grab users from database
    user0 = User.query.filter(User.name == user0_name).one()
    user1 = User.query.filter(User.name == user1_name).one()

    # grabbing vectors from each tweet in user.tweets
    user0_vects = np.array([tweet.vect for tweet in user0.tweets])
    user1_vects = np.array([tweet.vect for tweet in user1.tweets])

    # vertically stack (combine vectors) to train model
    vects = np.vstack([user0_vects, user1_vects])

    # generate labels for vects array
    labels = np.concatenate([np.zeros(len(user0.tweets)),
                             np.ones(len(user1.tweets))])

    # instantiate and train model
    log_reg = LogisticRegression().fit(vects, labels)

    # use nlp to generate embeddings - vectorize tweet()
    hypo_tweet_vect = vectorize_tweet(hypo_tweet_text).reshape(1, -1)

    # predicts and returns 0 or 1 depending upone logistic regression models prediction
    return log_reg.predict(hypo_tweet_vect)

max_stars_repo_name: BobBriksz/twitoff-bobbriksz

max_issues_repo_name: BobBriksz/twitoff-bobbriksz

max_forks_repo_name: BobBriksz/twitoff-bobbriksz


text:
# import the necessary packages
from libraries.videostream import WebcamVideoStream as vstream
import numpy as np
import argparse
import time
import cv2
import os
os.environ['TF_CPP_MIN_LOG_LEVEL'] = '3'
import sys
import time
import libraries.plate_detector as p_reader

# construct the argument parser and parse the arguments
original = sys.stdout
ap = argparse.ArgumentParser()
ap.add_argument("-m", "--model",default="tflite_model.tflite",
	help="path to tflite trained model [STRING]")
ap.add_argument("-a", "--address",default=0,
	help="stream address\n[rtsp] = rtsp://[username]:[password]@[ipaddress]/[camera specific additions] [STRING]\n[webcam] = 0 [INT]")
ap.add_argument("-cl", "--classes",default="obj.names",
	help="path to classes file of the model [STRING]")
ap.add_argument("-l", "--log",type=bool,default=True,
	help="creating a log of recognized numbers [True\False]")
ap.add_argument("-v", "--view",type=bool,default=False,
	help="viewing result for debugging [True/False]")
ap.add_argument("-sl", "--sleep",type=int,default=0,
	help="amount of sleep time between frames in seconds [DOUBLE]")
ap.add_argument("-si", "--size",type=int,default=416,
	help="size of the pictures the model was trained on [INT]")
ap.add_argument("-c", "--confidence", type=float, default=0.3,
	help="minimum probability to filter weak detections [DOUBLE]")
ap.add_argument("-sig", "--sigma",  default="ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789",
	help="the alphabet allowed in the licence plate [STRING]")
ap.add_argument("-lr", "--lengthrange", default="0",
	help="the length allowed for a licence plate {format: 7-8 or 3-6 etc.} [STRING]")
args = vars(ap.parse_args())

# loads tensorflow and the detector library
print("[INFO] loading detector...")
import libraries.tf_lite_detect as detector

# loads the classes file
print("[INFO] processing the classes file...")
classes = detector.read_class_names(args["classes"])

# initializing the plate reader
print("[INFO] loading plate detector...")
anpr = p_reader.plate_reader(debug=False)

print("[INFO] loading tesseract ocr...")
import pytesseract
pytesseract.pytesseract.tesseract_cmd = r'C:\Program Files\Tesseract-OCR\tesseract.exe'

# initialize the video stream, then allow the camera sensor to warm up
print("[INFO] starting video stream...")
if args["address"]=="0":
	args["assress"] = 0
vs = vstream(src=args["address"]).start()
time.sleep(1)
print("[INFO] Running...")
# start the FPS throughput estimator
fps = 0.0
startTime = time.time()
start = time.time()
end = time.time()

size=args["size"]
conf=args["confidence"]
model=args["model"]
view=args["view"]
sl= args["sleep"]
log= args["log"]
if args["lengthrange"]=="0":
	args["lengthrange"] = 0

log_file=open("log.txt","a")

# loop over frames from the video file stream
while (True):
	# grab the frame from the threaded video stream
	frame = vs.read()
	# construct a blob from the image
	s=time.time()
	boxes = detector.detect(framework='tflite', weights=model, size=size, tiny=True, model='yolov4',
							str_input=False, image=frame, pic_output=False, output='result.png', iou=0.45, score=conf,
							classes_file=classes)
	e=time.time()
	#print("tflite:",e-s)
	s=time.time()
	for x in boxes:
		if (x['class'] == "Car") or (x['class'] == "Motorcycle"):
			image = frame[x['c0']:x['c2'], x['c1']:x['c3']]
			(lpText, lpCnt) = anpr.find_and_ocr(image,args["sigma"],args["lengthrange"], psm=7,
												clearBorder=True)
			# only continue if the license plate was successfully OCR'd
			if lpText is not None and lpCnt is not None and lpText != "":
				print("[INFO] ", lpText)
				log_file.write(str(time.ctime())+"  --  "+lpText+"\n")
			if view and lpText is not None and lpCnt is not None and lpText != "":
				# fit a rotated bounding box to the license plate contour and
				# draw the bounding box on the license plate
				box = cv2.boxPoints(cv2.minAreaRect(lpCnt))
				box = box.astype("int")
				cv2.drawContours(image, [box], -1, (0, 0, 255), 2)
				# compute a normal (unrotated) bounding box for the license
				# plate and then draw the OCR'd license plate text on the
				# image
				(x, y, w, h) = cv2.boundingRect(lpCnt)
				cv2.putText(image, lpText, (x, y - 15),
							cv2.FONT_HERSHEY_SIMPLEX, 0.75, (0, 0, 255), 2)
				# show the output ANPR image
				cv2.imshow("Output ANPR", image)
				cv2.waitKey(0)

	e=time.time()
	#print("plates:",e-s)
	# update the FPS counter
	fps+=1.0
	# update the end time
	end = time.time()
	if(fps%100==0):
		print("[INFO] approx. FPS: {:.2f}".format(fps/(end-start)))
	key = cv2.waitKey(1) & 0xFF
	# if the `q` key was pressed, break from the loop
	if key == ord("q"):
		break
	time.sleep(sl)

# stop the timer and display FPS information
elapsed = time.time()-startTime
fps = fps/elapsed
print("[INFO] elasped time: {:.2f}".format(elapsed))
print("[INFO] approx. FPS: {:.2f}".format(fps))

# do a bit of cleanup
cv2.destroyAllWindows()
vs.stop()


text:
import numpy as np
import pytest

import reciprocalspaceship as rs

@pytest.mark.parametrize(
    "data",
    [
        np.linspace(0, 100, 1000),  # Simulates merged data
        np.repeat(np.linspace(0, 100, 1000), 3),  # Simulates unmerged data
    ],
)
@pytest.mark.parametrize("bins", [10, 20, 50])
@pytest.mark.parametrize("ascending", [True, False])
@pytest.mark.parametrize("shuffle", [True, False])
def test_bin_by_percentile(data, bins, ascending, shuffle):
    # Order shouldn't matter
    if shuffle:
        np.random.shuffle(data)

    assignments, labels = rs.utils.bin_by_percentile(data, bins, ascending)

    # Make sure number of bins is correct
    assert len(labels) == bins
    assert len(assignments) == len(data)
    assert len(np.unique(assignments)) == bins

    if ascending:
        np.all(np.diff(assignments) >= 0)
    else:
        np.all(np.diff(assignments) <= 0)

max_stars_repo_name: Hekstra-Lab/reciprocalspaceship

max_issues_repo_name: Hekstra-Lab/reciprocalspaceship

max_forks_repo_name: Hekstra-Lab/reciprocalspaceship


text:
# -*- coding: utf-8 -*-
"""
Created on Fri Dec 11 23:39:22 2015

@author: ajaver
"""

import os
import warnings

import numpy as np
import pandas as pd
import tables
from scipy.stats import chi2
from sklearn.covariance import MinCovDet

from tierpsy.helper.params import min_num_skel_defaults
from tierpsy.helper.misc import TimeCounter, print_flush, TABLE_FILTERS, save_modified_table, get_base_name

warnings.filterwarnings('ignore', '.*det > previous_det*',)
np.seterr(invalid='ignore')

worm_partitions = {'neck': (8, 16),
                   'midbody': (16, 33),
                   'hips': (33, 41),
                   # refinements of ['head']
                   'head_tip': (0, 4),
                   'head_base': (4, 8),    # ""
                   # refinements of ['tail']
                   'tail_base': (40, 45),
                   'tail_tip': (45, 49)}

name_width_fun = lambda part: 'width_' + part

def getValidIndexes(
        trajectories_data_f,
        min_num_skel=100,
        bad_seg_thresh=0.8,
        min_displacement=5,
        worm_index_type = 'worm_index_joined'):
    # min_num_skel - ignore trajectories that do not have at least this number of skeletons
    # min_dist - minimum distance explored by the blob to be consider a real
    # worm

    trajectories_data = trajectories_data_f[trajectories_data_f[worm_index_type] > 0]
    valid_ind_str = 'is_good_skel' if 'is_good_skel' in trajectories_data else 'has_skeleton'

    if len(trajectories_data[worm_index_type].unique()) == 1:
        good_skel_row = trajectories_data['skeleton_id'][
            trajectories_data[valid_ind_str].values.astype(np.bool)].values
        return (np.array([1]), good_skel_row)

    else:

        # get the fraction of worms that were skeletonized per trajectory
        how2agg = {
            valid_ind_str: [
                'mean', 'sum'], 'coord_x': [
                'max', 'min', 'count'], 'coord_y': [
                'max', 'min']}
        tracks_data = trajectories_data.groupby(
            worm_index_type).agg(how2agg)

        delX = tracks_data['coord_x']['max'] - tracks_data['coord_x']['min']
        delY = tracks_data['coord_y']['max'] - tracks_data['coord_y']['min']

        # /tracks_data['coord_x']['count']
        max_avg_dist = np.sqrt(delX * delX + delY * delY)

        skeleton_fracc = tracks_data[valid_ind_str]['mean']
        skeleton_tot = tracks_data[valid_ind_str]['sum']

        good_worm = (
            skeleton_fracc >= bad_seg_thresh) & (
            skeleton_tot >= min_num_skel)
        good_worm = good_worm & (max_avg_dist > min_displacement)

        good_traj_index = good_worm[good_worm].index

        good_row = (trajectories_data[worm_index_type].isin(good_traj_index)) \
            & (trajectories_data[valid_ind_str].values.astype(np.bool))
        good_skel_row = trajectories_data.loc[good_row, 'skeleton_id'].values

        return (good_traj_index, good_skel_row)

def _h_nodes2Array(skeletons_file, nodes4fit, valid_index=-1):
    '''
    Read the groups in skeletons file and save them as a matrix.
    Used by _h_readFeat2Check
    '''
    with tables.File(skeletons_file, 'r') as fid:
        assert all(node in fid for node in nodes4fit)

        if isinstance(valid_index, (float, int)) and valid_index < 0:
            valid_index = np.arange(fid.get_node(nodes4fit[0]).shape[0])

        n_samples = len(valid_index)
        n_features = len(nodes4fit)

        X = np.zeros((n_samples, n_features))

        if valid_index.size > 0:
            for ii, node in enumerate(nodes4fit):
                X[:, ii] = fid.get_node(node)[valid_index]

        return X

def _h_readFeat2Check(skeletons_file, valid_index=-1):
    nodes4fit = ['/skeleton_length', '/contour_area', '/width_midbody']

    worm_morph = _h_nodes2Array(skeletons_file, nodes4fit, valid_index)

    with tables.File(skeletons_file, 'r') as fid:
        if isinstance(valid_index, (float, int)) and valid_index < 0:
            valid_index = np.arange(fid.get_node(nodes4fit[0]).shape[0])

        cnt_widths = fid.get_node('/contour_width')[valid_index, :]

        sample_N = cnt_widths.shape[1]
        ht_limits = int(round(sample_N / 7))

        head_widths = cnt_widths[:, 1:ht_limits]
        tail_widths = cnt_widths[:, -ht_limits:-1]

        # The log(x+1e-1) transform skew the distribution to the right, so the lower values have a higher change to
        # be outliers.I do this because a with close to zero is typically an
        # oulier.
        head_widthsL = np.log(head_widths + 1)
        tail_widthsL = np.log(tail_widths + 1)

    return worm_morph, head_widths, tail_widths, head_widthsL, tail_widthsL

def _h_getMahalanobisRobust(dat, critical_alpha=0.01, good_rows=np.zeros(0)):
    '''Calculate the Mahalanobis distance from the sample vector.'''
    if good_rows.size == 0:
        good_rows = np.any(~np.isnan(dat), axis=1)

    try:
        dat2fit = dat[good_rows]
        assert not np.any(np.isnan(dat2fit))

        robust_cov = MinCovDet().fit(dat2fit)
        mahalanobis_dist = np.sqrt(robust_cov.mahalanobis(dat))
    except ValueError:
        # this step will fail if the covariance matrix is not singular. This happens if the data is not
        # a unimodal symetric distribution. For example there is too many small noisy particles. Therefore
        # I will take a safe option and return zeros in the mahalanobis
        # distance if this is the case.
        mahalanobis_dist = np.zeros(dat.shape[0])

    # critial distance of the maholanobis distance using the chi-square distirbution
    # https://en.wikiversity.org/wiki/Mahalanobis%27_distance
    # http://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.chi2.html
    maha_lim = chi2.ppf(1 - critical_alpha, dat.shape[1])
    outliers = mahalanobis_dist > maha_lim

    return mahalanobis_dist, outliers, maha_lim

def _h_getPerpContourInd(
        skeleton,
        skel_ind,
        contour_side1,
        contour_side2,
        contour_width):
    # get the closest point in the contour from a line perpedicular to the skeleton.
    #%%

    # get the slop of a line perpendicular to the keleton
    dR = skeleton[skel_ind + 1] - skeleton[skel_ind - 1]
    #m = dR[1]/dR[0]; M = -1/m
    a = -dR[0]
    b = +dR[1]

    c = b * skeleton[skel_ind, 1] - a * skeleton[skel_ind, 0]

    max_width_squared = np.max(contour_width)**2
    # modified from https://en.wikipedia.org/wiki/Distance_from_a_point_to_a_line
    #a = M, b = -1

    # make sure we are not selecting a point that get traversed by a coiled
    # worm
    dist2cnt1 = np.sum((contour_side1 - skeleton[skel_ind])**2, axis=1)
    d1 = np.abs(a * contour_side1[:, 0] - b * contour_side1[:, 1] + c)
    d1[dist2cnt1 > max_width_squared] = np.nan

    dist2cnt2 = np.sum((contour_side2 - skeleton[skel_ind])**2, axis=1)
    d2 = np.abs(a * contour_side2[:, 0] - b * contour_side2[:, 1] + c)
    d2[dist2cnt2 > max_width_squared] = np.nan

    try:
        cnt1_ind = np.nanargmin(d1)
        cnt2_ind = np.nanargmin(d2)
    except ValueError:
        cnt1_ind = np.nan
        cnt2_ind = np.nan

    return cnt1_ind, cnt2_ind

def _h_calcArea(cnt):
    # make sure the contours are in the counter-clockwise direction
    # x1y2 - x2y1(http://mathworld.wolfram.com/PolygonArea.html)
    signed_area = np.sum(cnt[:-1, 0] * cnt[1:, 1] - cnt[1:, 0] * cnt[:-1, 1])
    return np.abs(signed_area / 2)

def filterPossibleCoils(
        skeletons_file,
        max_width_ratio=2.25,
        max_area_ratio=6):
    with pd.HDFStore(skeletons_file, 'r') as table_fid:
        trajectories_data = table_fid['/trajectories_data']

    if not 'is_good_skel' in trajectories_data:
        is_good_skel = trajectories_data['has_skeleton'].values.copy()
    else:
        is_good_skel = trajectories_data['is_good_skel'].values.copy()

    tot_skeletons = len(is_good_skel)
    with tables.File(skeletons_file, 'r') as ske_file_id:

        skeletons = ske_file_id.get_node('/skeleton')
        contour_side1s = ske_file_id.get_node('/contour_side1')
        contour_side2s = ske_file_id.get_node('/contour_side2')
        contour_widths = ske_file_id.get_node('/contour_width')

        sample_N = contour_widths.shape[1]

        ht_limits = int(round(sample_N / 6))
        mid_limits = (int(round(3 * sample_N / 6)),
                      int(round(4 * sample_N / 6)) + 1)

        for skel_id in range(tot_skeletons):
            if is_good_skel[skel_id] == 0:
                continue

            contour_width = contour_widths[skel_id]
            contour_side1 = contour_side1s[skel_id]
            contour_side2 = contour_side2s[skel_id]
            skeleton = skeletons[skel_id]

            #%%
            edge_length = 8
            p1 = contour_side1[:-edge_length]
            p2 = contour_side1[edge_length:]

            points = contour_side1[edge_length // 2:-edge_length // 2]
            ang2 = np.arctan2(points[:, 0] - p2[:, 0], points[:, 1] - p2[:, 1])
            ang1 = np.arctan2(p1[:, 0] - points[:, 0], p1[:, 1] - points[:, 1])
            angles = ang2 - ang1

            for i in range(angles.size):
                if angles[i] > np.pi:
                    angles[i] = angles[i] - 2 * np.pi
                elif angles[i] < -np.pi:
                    angles[i] = angles[i] + 2 * np.pi
            angles = angles * 180 / np.pi

            blurWin = np.full(edge_length, 1. / edge_length)
            anglesb = np.convolve(angles, blurWin, 'same')
            #%%
            # the idea is that when the worm coils and there is an skeletons, it is
            # likely to be a cnonsequence of the head/tail protuding, therefore we can
            # use the head/tail withd to get a good ratio of the worm width

            # calculate head and tail width
            head_w = contour_width[ht_limits]
            tail_w = contour_width[-ht_limits]
            midbody_w = np.max(contour_width)

            '''
            Does the worm more than double its width from the head/tail?
            Note: if the worm coils, its width will grow to more than
            double that at the end of the head.
            '''
            if midbody_w / head_w > max_width_ratio or midbody_w / \
                    tail_w > max_width_ratio or max(head_w / tail_w, tail_w / head_w) > max_width_ratio:
                is_good_skel[skel_id] = 0
                continue

            # calculate the head and tail area (it is an approximation the
            # limits are not super well defined, but it is enough for
            # filtering)
            head_skel_lim = skeleton[np.newaxis, ht_limits]
            tail_skel_lim = skeleton[np.newaxis, -ht_limits]

            cnt_side1_ind_h, cnt_side2_ind_h = _h_getPerpContourInd(
                skeleton, ht_limits, contour_side1, contour_side2, contour_width)
            cnt_side1_ind_t, cnt_side2_ind_t = _h_getPerpContourInd(
                skeleton, -ht_limits, contour_side1, contour_side2, contour_width)

            if cnt_side1_ind_h > cnt_side1_ind_t or cnt_side2_ind_h > cnt_side2_ind_t or \
            np.any(np.isnan([cnt_side1_ind_h, cnt_side2_ind_h, cnt_side1_ind_t, cnt_side2_ind_t])):
                is_good_skel[skel_id] = 0
                continue

            # if cnt_side1_ind_h>cnt_side1_ind_t:
            #    cnt_side1_ind_h, cnt_side1_ind_t = cnt_side1_ind_t, cnt_side1_ind_h
            # if cnt_side2_ind_h>cnt_side2_ind_t:
            #    cnt_side2_ind_h, cnt_side2_ind_t = cnt_side2_ind_t, cnt_side2_ind_h

            cnt_head = np.concatenate((contour_side1[:cnt_side1_ind_h + 1], head_skel_lim,
                                           contour_side2[:cnt_side2_ind_h + 1][::-1]))

            cnt_tail = np.concatenate(
                (contour_side2[cnt_side2_ind_t:][::-1],
                    tail_skel_lim,
                    contour_side1[cnt_side1_ind_t:]))

            area_head = _h_calcArea(cnt_head)
            area_tail = _h_calcArea(cnt_tail)

            '''Is the tail too small (or the head too large)?
            Note: the area of the head and tail should be roughly the same size.
            A 2-fold difference is huge!
            '''
            if area_tail == 0 or area_head == 0 or area_head / \
                    area_tail > max_width_ratio or area_tail / area_head > max_width_ratio:
                is_good_skel[skel_id] = 0
                continue

            # calculate the area of the rest of the body
            cnt_rest = np.concatenate(
                (head_skel_lim,
                 contour_side1[
                     cnt_side1_ind_h:cnt_side1_ind_t + 1],
                    tail_skel_lim,
                    contour_side2[
                     cnt_side2_ind_h:cnt_side2_ind_t + 1][
                     ::-1],
                    head_skel_lim))
            area_rest = _h_calcArea(cnt_rest)

            '''
            Are the head and tail too small (or the body too large)?
            Note: earlier, the head and tail were each chosen to be 4/24 = 1/6
            the body length of the worm. The head and tail are roughly shaped
            like rounded triangles with a convex taper. And, the width at their
            ends is nearly the width at the center of the worm. Imagine they were
            2 triangles that, when combined, formed a rectangle similar to the
            midsection of the worm. The area of this rectangle would be greater
            than a 1/6 length portion from the midsection of the worm (the
            maximum area per length in a worm is located at its midsection). The
            combined area of the right and left sides is 4/6 of the worm.
            Therefore, the combined area of the head and tail must be greater
            than (1/6) / (4/6) = 1/4 the combined area of the left and right
            sides.
            '''
            if area_rest / (area_head + area_tail) > max_area_ratio:
                is_good_skel[skel_id] = 0
                continue
    trajectories_data['is_good_skel'] = is_good_skel.astype(np.uint8)
    save_modified_table(skeletons_file, trajectories_data, 'trajectories_data')

def _h_calAreaSignedArray(cnt_side1, cnt_side2):
    '''calculate the contour area using the shoelace method, the sign indicate the contour orientation.'''
    assert cnt_side1.shape == cnt_side2.shape
    if cnt_side1.ndim == 2:
        # if it is only two dimenssion (as if in a single skeleton).
        # Add an extra dimension to be compatible with the rest of the code
        cnt_side1 = cnt_side1[np.newaxis, :, :]
        cnt_side2 = cnt_side2[np.newaxis, :, :]

    contour = np.hstack((cnt_side1, cnt_side2[:, ::-1, :]))
    signed_area = np.sum(
        contour[:,:-1,0] * contour[:,1:,1] -
        contour[:,1:,0] * contour[:,:-1,1],
        axis=1)/ 2

    assert signed_area.size == contour.shape[0]
    return signed_area

def _h_calAreaArray(cnt_side1, cnt_side2):
    '''calculate the contour area using the shoelace method'''
    signed_area = _h_calAreaSignedArray(cnt_side1, cnt_side2)
    return np.abs(signed_area)

def _addMissingFields(skeletons_file):
    '''
    Add missing fields that might have not been calculated before in older videos.
    '''
    with tables.File(skeletons_file, 'r+') as fid:
        if not '/contour_area' in fid:
            contour_side1 = fid.get_node('/contour_side1')[:]
            contour_side2 = fid.get_node('/contour_side2')[:]
            cnt_area = _h_calAreaArray(contour_side1, contour_side2)
            fid.create_carray('/', "contour_area", obj=cnt_area, filters=TABLE_FILTERS)

        if not '/width_midbody' in fid:
            midbody_ind = (17, 32)
            contour_width = fid.get_node('/contour_width')[:]
            width_midbody = np.median(contour_width[:, midbody_ind[0]:midbody_ind[1]+1], axis=1)
            fid.create_carray('/', "width_midbody", obj=width_midbody, filters=TABLE_FILTERS)

def filterByPopulationMorphology(skeletons_file, good_skel_row, critical_alpha=0.01):
    base_name = get_base_name(skeletons_file)
    progress_timer = TimeCounter('')

    print_flush(base_name + ' Filter Skeletons: Starting...')
    with pd.HDFStore(skeletons_file, 'r') as table_fid:
        trajectories_data = table_fid['/trajectories_data']

    if not 'is_good_skel' in trajectories_data:
        trajectories_data['is_good_skel'] = trajectories_data['has_skeleton']

    if good_skel_row.size > 0:
        # nothing to do if there are not valid skeletons left.

        print_flush(
            base_name +
            ' Filter Skeletons: Reading features for outlier identification.')

        #add possible missing fields that were con calculated in older versions of the software
        _addMissingFields(skeletons_file)

        # calculate classifier for the outliers
        nodes4fit = ['/skeleton_length', '/contour_area', '/width_midbody']
        worm_morph = _h_nodes2Array(skeletons_file, nodes4fit, -1)
        #worm_morph[~trajectories_data['is_good_skel'].values] = np.nan
        feats4fit = [worm_morph]

        #feats4fit = _h_readFeat2Check(skeletons_file)

        print_flush(
            base_name +
            ' Filter Skeletons: Calculating outliers. Total time:' +
            progress_timer.get_time_str())

        tot_rows2fit = feats4fit[0].shape[0]
        # check all the data to fit has the same size in the first axis
        assert all(tot_rows2fit == featdat.shape[0] for featdat in feats4fit)
        outliers_rob = np.zeros(tot_rows2fit, np.bool)
        outliers_flag = np.zeros(tot_rows2fit, np.int)
        assert len(feats4fit) < 64  # otherwise the outlier flag will not work

        for out_ind, dat in enumerate(feats4fit):
            maha, out_d, lim_d = _h_getMahalanobisRobust(
                dat, critical_alpha, good_skel_row)
            outliers_rob = outliers_rob | out_d

            # flag the outlier flag by turning on the corresponding bit
            outliers_flag += (out_d) * (2**out_ind)

        print_flush(
            base_name +
            ' Filter Skeletons: Labeling valid skeletons. Total time:' +
            progress_timer.get_time_str())

        # labeled rows of valid individual skeletons as GOOD_SKE
        trajectories_data['is_good_skel'] &= ~outliers_rob
        trajectories_data['skel_outliers_flag'] = outliers_flag

    # Save the new is_good_skel column
    if trajectories_data['is_good_skel'].dtypes == bool:
        trajectories_data['is_good_skel'] = trajectories_data[
            'is_good_skel'].astype(np.uint8)
    save_modified_table(skeletons_file, trajectories_data, 'trajectories_data')

    print_flush(
        base_name +
        ' Filter Skeletons: Finished. Total time:' +
        progress_timer.get_time_str())

def getFilteredSkels(
        skeletons_file,
        min_num_skel=100,
        bad_seg_thresh=0.8,
        min_displacement=5,
        critical_alpha=0.01,
        max_width_ratio=2.25,
        max_area_ratio=6):

    min_num_skel = min_num_skel_defaults(skeletons_file, min_num_skel=min_num_skel)

    # check if the skeletonization finished succesfully
    with tables.File(skeletons_file, "r") as ske_file_id:
        skeleton_table = ske_file_id.get_node('/skeleton')

    #eliminate skeletons that do not match a decent head, tail and body ratio. Likely to be coils. Taken from Segworm.
    filterPossibleCoils(
        skeletons_file,
        max_width_ratio=max_width_ratio,
        max_area_ratio=max_area_ratio)

    with pd.HDFStore(skeletons_file, 'r') as table_fid:
        trajectories_data = table_fid['/trajectories_data']

    # get valid rows using the trajectory displacement and the
    # skeletonization success. These indexes will be used to calculate statistics of what represent a valid skeleton.
    good_traj_index, good_skel_row = getValidIndexes(
        trajectories_data, min_num_skel=min_num_skel, bad_seg_thresh=bad_seg_thresh, min_displacement=min_displacement)

    #filter skeletons depending the population morphology (area, width and length)
    filterByPopulationMorphology(
        skeletons_file,
        good_skel_row,
        critical_alpha=critical_alpha)

max_stars_repo_name: JGNieto/tierpsy-tracker

max_issues_repo_name: RongkangXiong/tierpsy-tracker

max_forks_repo_name: RongkangXiong/tierpsy-tracker


text:
import torch
import torch.autograd
from torch.autograd import Variable
import torch.optim as optim
import torch.nn as nn

import numpy as np
import scipy as sp
import scipy.linalg

from qpsolvers import solve_qp

from core.DDPG import DDPGagent
from core.ConstraintNetwork import ConstraintNetwork

class SafeDDPGagent(DDPGagent):
    def __init__(self, state_dim, act_dim, constraint_dim, num_agents, constraint_networks_dir, col_margin = 0.35 ,
            hidden_size=256, actor_learning_rate=1e-4, critic_learning_rate=1e-3,
            gamma=0.99, tau=1e-2, max_memory_size=16000,  soften = True):

        # Call DDPGagent's constructor
        super().__init__(state_dim, act_dim, num_agents ,hidden_size,
                         actor_learning_rate, critic_learning_rate, gamma,
                         tau, max_memory_size)

        # Extra Params
        self.col_margin = col_margin
        self.constraint_dim = constraint_dim

        self.total_state_dim = self.state_dim * self.num_agents
        self.total_constraint_dim = self.constraint_dim * self.num_agents
        self.total_action_dim = self.act_dim * self.num_agents

        self.constraint_nets = self.total_constraint_dim*[None]

        # Initialize constraint networks
        for i in range(self.total_constraint_dim):
            self.constraint_nets[i] = ConstraintNetwork(self.total_state_dim, self.total_action_dim)
            self.constraint_nets[i].load_state_dict(torch.load(constraint_networks_dir
                                                    + "constraint_net_" + str(i) + ".pkl"))

        # Define Solver Globaly
        self.solver_interventions = 0
        self.solver_infeasible    = 0

        # Choose Solver
        if soften:
            self.correct_actions = self.correct_actions_soften
        else:
            self.correct_actions = self.correct_actions_hard

    def reset_metrics(self):
        self.solver_interventions = 0
        self.solver_infeasible    = 0

    def get_interventions(self):
        return self.solver_interventions

    def get_infeasible(self):
        return self.solver_infeasible

    @torch.no_grad()
    def get_action(self, state, constraint):
        state = Variable(torch.from_numpy(state).float().unsqueeze(0))
        action = self.actor(state)
        action = action.detach().numpy().flatten()

        # transform numpy array into list of 3 actions
        action = self.correct_actions(state, action, constraint)

        actions = np.split(action, self.num_agents)
        return actions

    @torch.no_grad()
    def correct_actions_hard(self, state, actions, constraint):

        # (1) Problem Variables
        # Problem specific constants
        I    = np.eye(self.total_action_dim)
        ones = np.ones(self.total_action_dim)
        C    = np.concatenate(constraint)

        # Formulate the constraints using neural networks
        G    = np.zeros([self.total_action_dim, self.total_action_dim])
        for i, net in enumerate(self.constraint_nets):
            G[i, :] = net(state).numpy()

        # (2) Problem Variables in QP form
        # Cost Function
        q = -actions
        P = np.eye(self.total_action_dim)

        # Constraints
        A = np.concatenate([-G, I, -I])
        ub = np.concatenate([C - self.col_margin, ones, ones])
        lb = None

        # Solve Optimization Problem
        try:
            x = solve_qp(P.astype(np.float64), q.astype(np.float64), A.astype(np.float64),
                         ub.astype(np.float64), None, None, None, None)
        except:
            self.solver_infeasible +=1
            return actions

        # Count Solver interventions
        if np.linalg.norm(actions - x) > 1e-3:
            self.solver_interventions += 1

        return x

    @torch.no_grad()
    def correct_actions_soften(self, state, actions, constraint):

        # (1) Create solver as a globar variable
        l1_penalty = 1000

        # (2) Problem Variables
        # Problem specific constants
        I     = np.eye(self.total_action_dim)
        Z     = np.zeros([self.total_action_dim, self.total_action_dim])
        ones  = np.ones(self.total_action_dim)
        zeros = np.zeros(self.total_action_dim)
        C     = np.concatenate(constraint) - self.col_margin

        # Formulate the constraints using neural networks
        G    = np.zeros([self.total_action_dim, self.total_action_dim])
        for i, net in enumerate(self.constraint_nets):
            G[i, :] = net(state).numpy()

        # (2) Problem Variables in QP form
        # Cost Function
        P = sp.linalg.block_diag(I, Z + I * 0.001, Z + I * 0.001)
        q = np.concatenate([-actions, ones, zeros])

        # Constraints
        A = np.vstack((np.concatenate([-G, Z, -I], axis = 1),
                       np.concatenate([Z, Z, -I], axis = 1),
                       np.concatenate([Z, -I,  l1_penalty * I], axis = 1),
                       np.concatenate([Z, -I, -l1_penalty * I], axis = 1)))

        ub = np.concatenate((C, zeros, zeros, zeros))
        lb = None

        # Solve Optimization Problem
        try:
            x = solve_qp(P.astype(np.float64), q.astype(np.float64), A.astype(np.float64),
                         ub.astype(np.float64), None, None, None, None)

            x = x[0:(self.total_action_dim)]
        except:
            self.solver_infeasible +=1
            return actions

        # Count Solver interventions
        if np.linalg.norm(actions - x) > 1e-3:
            self.solver_interventions += 1

        return x

max_stars_repo_name: zisikons/deep-rl

max_issues_repo_name: zisikons/deep-rl

max_forks_repo_name: zisikons/deep-rl

reason: import numpy,import scipy


text:
import matplotlib.pyplot as plt
import numpy as np
import csv
from sklearn import metrics
import sys

log_loc0 = 'logs/comparisons/dnn/dim-comparison/'
log_loc1 = 'logs/comparisons/qgnn/iteration_comparison/'
log_loc2 = 'logs/comparisons/qgnn/learning_rate_comparison/'
log_loc3 = 'logs/comparisons/qgnn/dimension_comparison/'
log_loc4 = 'logs/comparisons/qgnn/arch_comparison/'

run_type = 'test'

# TO DO:
# some run types are deprecated, need to update them

# TO DO:
# implement exectuable input form
# ask run_type as input

if run_type == 'cla-qua':
	log_location0 = log_loc0 + 'hid1/'
	log_location1 = log_loc0 + 'hid10/'
	log_location2 = log_loc0 + 'hid100/'
	log_location3 = log_loc1 + 'it1/'

	pdf_location = 'pdf/comparison/cla-qua/'
	png_location = 'png/comparison/cla-qua/'

	label0 = r'$Classical \, Network \, with \, N_{Dim} = 1$'
	label1 = r'$Classical \, Network \, with \, N_{Dim} = 10$'
	label2 = r'$Classical \, Network \, with \, N_{Dim} = 100$'
	label3 = r'$Quantum \, Network \, with \, N_{Dim} = 1$'

if run_type == 'classical_dimension_comparison':
	log_location0 = log_loc + 'hid1/'
	log_location1 = log_loc + 'hid10/'
	log_location2 = log_loc + 'hid100/'

	label1 = r'$N_{dim}=1$'
	label2 = r'$N_{dim}=10$'
	label3 = r'$N_{dim}=100$'

if run_type == 'iteration_comparison':
	# Iteration comparison
	n_runs = 4
	log_location0 = log_loc1 + 'it1/'
	log_location1 = log_loc1 + 'it2/'
	log_location2 = log_loc1 + 'it3/'
	log_location3 = log_loc1 + 'it4/'

	pdf_location = 'pdf/comparison/qgnn/iteration_comparison/'
	png_location = 'png/comparison/qgnn/iteration_comparison/'

	label0 = r'$N_{it}=1$'
	label1 = r'$N_{it}=2$'
	label2 = r'$N_{it}=3$'
	label3 = r'$N_{it}=4$'

if run_type == 'l.rate_comparison':
	# Learning rate comparison
	n_runs = 3
	log_location0 = log_loc2 + 'lr_1e-2/'
	log_location1 = log_loc2 + 'lr_2e-2/'
	log_location2 = log_loc2 + 'lr_3e-2/'

	pdf_location = 'pdf/comparison/qgnn/learning_rate_comparison/'
	png_location = 'png/comparison/qgnn/learning_rate_comparison/'

	label0 = r'$lr=1x10^{-2}$'
	label1 = r'$lr=2x10^{-2}$'
	label2 = r'$lr=3x10^{-2}$'

if run_type == 'dimension_comparison':
	## Dimension comparison
	n_runs = 2
	log_location0 = log_loc2 + 'lr_3e-2/'
	log_location1 = log_loc3 + 'dim2/'

	pdf_location = 'pdf/comparison/qgnn/dimension_comparison/'
	png_location = 'png/comparison/qgnn/dimension_comparison/'

	label0 = r'$N_{dim}=1$'
	label1 = r'$N_{dim}=2$'

if run_type == 'test':## Test comparison
	n_runs = 2
	log_location0 = log_loc2 + 'lr_3e-2/'
	log_location1 = 'logs/test/'

	pdf_location = 'pdf/test/'
	png_location = 'png/test/'

	label0 = 'base'
	label1 = 'test'

if run_type == 'arch_comparison':
	## General vs. Real Rotations
	n_runs = 2
	log_location0 = log_loc2 + 'lr_3e-2/'
	log_location1 = log_loc4 + 'embedding_learning/'

	pdf_location = 'pdf/comparison/qgnn/arch_comparison/'
	png_location = 'png/comparison/qgnn/arch_comparison/'

	label0 = 'Old Architecture'
	label1 = 'New(test) Architecture'

print('Comparison plots will be saved to:')
print('PDF: ' + pdf_location)
print('PNG: ' + png_location)
print('Comparing %d runs................'%n_runs)

interval = 50

with open(log_location0+'log_validation.csv', 'r') as f:
	reader = csv.reader(f, delimiter=',')
	valid0 = np.array(list(reader)).astype(float)
	valid0 = valid0[:28]
x0 = [i*interval for i  in range(len(valid0))]
with open(log_location1+'log_validation.csv', 'r') as f:
	reader = csv.reader(f, delimiter=',')
	valid1 = np.array(list(reader)).astype(float)
	valid1 = valid1[:28]
x1 = [i*interval for i  in range(len(valid1))]
if n_runs > 2:
	with open(log_location2+'log_validation.csv', 'r') as f:
		reader = csv.reader(f, delimiter=',')
		valid2 = np.array(list(reader)).astype(float)
		valid2 = valid2[:28]
	x2 = [i*interval for i  in range(len(valid2))]
if n_runs > 3 :
	with open(log_location3+'log_validation.csv', 'r') as f:
		reader = csv.reader(f, delimiter=',')
		valid3 = np.array(list(reader)).astype(float)
		valid3 = valid3[:28]
	x3 = [i*interval for i  in range(len(valid3))]

plt.clf()
plt.plot(x0,valid0[:,0],label=label0,c='darkorange')
plt.plot(x1,valid1[:,0],label=label1,c='navy')
if n_runs > 2: plt.plot(x2,valid2[:,0],label=label2,c='red')
if n_runs > 3: plt.plot(x3,valid3[:,0],label=label3,c='green')
plt.title('Validation Accuracy')
plt.xlabel('Update')
plt.ylabel('Accuracy')
plt.legend()
plt.savefig(pdf_location+'validation_accuracy.pdf')
plt.savefig(png_location+'validation_accuracy.png')

plt.clf()
plt.plot(x0,valid0[:,1],label=label0,c='darkorange')
plt.plot(x1,valid1[:,1],label=label1,c='navy')
if n_runs > 2: plt.plot(x2,valid2[:,1],label=label2,c='red')
if n_runs > 3: plt.plot(x3,valid3[:,1],label=label3,c='green')
plt.title('Validation AUC')
plt.xlabel('Update')
plt.ylabel('AUC')
plt.legend()
plt.savefig(pdf_location+'validation_auc.pdf')
plt.savefig(png_location+'validation_auc.png')

plt.clf()
plt.plot(x0,valid0[:,2],label=label0,c='darkorange')
plt.plot(x1,valid1[:,2],label=label1,c='navy')
if n_runs > 2: plt.plot(x2,valid2[:,2],label=label2,c='red')
if n_runs > 3: plt.plot(x3,valid3[:,2],label=label3,c='green')
plt.title('Validation Loss')
plt.xlabel('Update')
plt.ylabel('Loss')
plt.legend()
plt.savefig(pdf_location+'validation_loss.pdf')
plt.savefig(png_location+'validation_loss.png')

with open(log_location0 + 'log_loss.csv', 'r') as f:
	reader = csv.reader(f, delimiter=',')
	loss0 = np.array(list(reader)).astype(float)
	loss0 = loss0[:1399]
x0 = [i+1 for i  in range(len(loss0))]

with open(log_location1 + 'log_loss.csv', 'r') as f:
	reader = csv.reader(f, delimiter=',')
	loss1 = np.array(list(reader)).astype(float)
	loss1 = loss1[:1399]
x1 = [i+1 for i  in range(len(loss1))]

if n_runs > 2:
	with open(log_location2 + 'log_loss.csv', 'r') as f:
		reader = csv.reader(f, delimiter=',')
		loss2 = np.array(list(reader)).astype(float)
		loss2 = loss2[:1399]
	x2 = [i+1 for i  in range(len(loss2))]

if n_runs > 3:
	with open(log_location3 + 'log_loss.csv', 'r') as f:
		reader = csv.reader(f, delimiter=',')
		loss3 = np.array(list(reader)).astype(float)
		loss3 = loss3[:1399]
	x3 = [i+1 for i  in range(len(loss3))]

plt.clf()
plt.plot(x0,loss0,label=label0,c='darkorange')
plt.plot(x1,loss1,label=label1,c='navy')
if n_runs > 2: plt.plot(x2,loss2,label=label2,c='red')
if n_runs > 3: plt.plot(x3,loss3,label=label3,c='green')
plt.title('Training Loss')
plt.xlabel('Update')
plt.ylabel('Loss')
plt.legend()
plt.savefig(pdf_location+'training_loss.pdf')
plt.savefig(png_location+'training_loss.png')

max_stars_repo_name: QC-SciLab/HepTrkX-quantum

max_issues_repo_name: QC-SciLab/HepTrkX-quantum

max_forks_repo_name: QC-SciLab/HepTrkX-quantum


text:
import codecs
import networkx as nx
from itertools import combinations

f = codecs.open('data/corpus/docs_topic.tsv', 'r', 'utf-8')
edges = []
nodes = [0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,
         25,26]
# name = {
#     '0': 'Nagyobb csoportok',
#     '1': 'Kisebb csoportok: egyéni portrék, stb.',
#     '2': 'Úton',
#     '3': 'Kisebb csoportok: családok, egyének',
#     '4': 'Egyéni portrék, politikusok, stb.',
#     '5': 'Adatok',
#     '6': 'Átkelés'
# }

name = {}
for i in range(0,27):
    name[str(i).zfill(2)] = 'Topic ' + str(i).zfill(2)

nodes = [str(n).zfill(2) for n in nodes]
edge_counts = {}
possible_edges = list(combinations(nodes, 2))
for edge in possible_edges:
    edge_counts[edge] = 0

for l in f:
    l = l.strip().split('\t')
    ns = sorted(l[1].split('|'))
    combos = list(combinations(ns, 2))
    for combo in combos:
        if combo in edge_counts.keys():
            edge_counts[combo] += 1
        elif (combo[1], combo[0]) in edge_counts.keys():
            edge_counts[(combo[1], combo[0])]

G = nx.Graph()
for n in nodes:
    G.add_node(n, name=name[n])

for e in edge_counts:
    G.add_edge(e[0], e[1], weight=edge_counts[e])

nx.write_graphml(G, "data/corpus/topic_graph.graphml")

max_stars_repo_name: k-monitor/parldata

max_issues_repo_name: k-monitor/parldata

max_forks_repo_name: k-monitor/parldata

reason: import networkx


text:
import argparse
import os

import IPython
import matplotlib.pyplot as plt
import numpy as np
import scipy as sp

import torch
from context import utils
from utils.misc import get_equal_dicts, length_of_longest
import utils.plotting as plot
from data_analysis import invert_signs

"""Script for analyzing experiment E005

This experiment examines the scaling of the parallel implementation by running
different numbers of perturbations on different numbers of CPUs.

From the submission script:

    PERTURBATIONS=(2 4 8 16 32 64 128 256 512 1024 2048 4096 8192 16384)
    CORES=(1 2 4 8 12 16 20 24)

Second run of jobs starts at job id 792529
"""

parser = argparse.ArgumentParser(description='Monitorer')
parser.add_argument('-d', type=str, default=None, metavar='--directory', help='The directory of checkpoints to monitor.')
args = parser.parse_args()

if not args.d:
    # args.d = "/home/jakob/mnt/Documents/es-rl/experiments/checkpoints/E001"
    args.d = "/home/jakob/Dropbox/es-rl/experiments/checkpoints/E001-SM"
    # args.d = "/Users/Jakob/mnt/Documents/es-rl/experiments/checkpoints/E001"

save_dir = os.path.join(args.d, 'analysis')
if not os.path.exists(save_dir):
    os.mkdir(save_dir)

directories = [os.path.join(args.d, di) for di in os.listdir(args.d) if os.path.isdir(os.path.join(args.d, di)) and di != 'monitoring']
filename = 'state-dict-algorithm.pkl'
rm_filenames = ['state-dict-best-algorithm.pkl', 'state-dict-best-optimizer.pkl', 'state-dict-best-model.pkl']

algorithm_states = []
workers = []
perturbations = []
for i, d in enumerate(directories):
    try:
        s = torch.load(os.path.join(d, filename))
        algorithm_states.append(s)
        print("Loaded " + str(i) + "/" + str(len(directories)) + ": " + d)
    except:
        print("No files found in: " + d)
invert_signs(algorithm_states, keys='all')

groups = [0]*len(algorithm_states)
for i, s in enumerate(algorithm_states):
    if s['safe_mutation'] == 'SUM':
        groups[i] = 1
groups = np.array(groups)

for k in ['return_unp', 'return_avg', 'return_min', 'return_max']:

    list_of_series = [s['stats'][k] for s in algorithm_states]
    list_of_genera = [s['stats']['generations'] for s in algorithm_states]
    l = length_of_longest(list_of_series)
    indices = [i for i, series in enumerate(list_of_series) if len(series) == l]

    list_of_series = [list_of_series[i] for i in indices]
    list_of_genera = [list_of_genera[i] for i in indices]
    groups_longest_series = groups[indices]

    plot.timeseries_median_grouped(list_of_genera, list_of_series, groups_longest_series, xlabel='generations', ylabel=k)
    plt.savefig(os.path.join(save_dir, 'analysis-01-' + k + '.pdf'), bbox_inches='tight')
    plt.close()

IPython.embed()

max_stars_repo_name: JakobHavtorn/es-rl

max_issues_repo_name: JakobHavtorn/es-rl

max_forks_repo_name: JakobHavtorn/es-rl