text:
# Copyright 2017 Battelle Energy Alliance, LLC
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
"""
  This module contains the Simultaneous Perturbation Stochastic Approximation Optimization strategy

  Created on June 16, 2016
  @ author: chenj, alfoa
"""
#for future compatibility with Python 3--------------------------------------------------------------
from __future__ import division, print_function, unicode_literals, absolute_import
import warnings
warnings.simplefilter('default',DeprecationWarning)
#End compatibility block for Python 3----------------------------------------------------------------

#External Modules------------------------------------------------------------------------------------
import sys
import os
import copy
import numpy as np
import scipy
#External Modules End--------------------------------------------------------------------------------

#Internal Modules------------------------------------------------------------------------------------
from .GradientBasedOptimizer import GradientBasedOptimizer
import Distributions
from utils import mathUtils,randomUtils,InputData
import SupervisedLearning
#Internal Modules End--------------------------------------------------------------------------------

class SPSA(GradientBasedOptimizer):
  """
    Simultaneous Perturbation Stochastic Approximation Optimizer
  """
  @classmethod
  def getInputSpecification(cls):
    """
      Method to get a reference to a class that specifies the input data for class cls.
      @ In, cls, the class for which we are retrieving the specification
      @ Out, inputSpecification, InputData.ParameterInput, class to use for specifying input of cls
    """
    inputSpecification = super(SPSA,cls).getInputSpecification()

    # add additional parameters to "parameter"
    param = inputSpecification.popSub('parameter')
    param.addSub(InputData.parameterInputFactory('gamma', contentType=InputData.FloatType, strictMode=True))
    param.addSub(InputData.parameterInputFactory('c'    , contentType=InputData.FloatType, strictMode=True))
    param.addSub(InputData.parameterInputFactory('a'    , contentType=InputData.FloatType, strictMode=True))
    param.addSub(InputData.parameterInputFactory('alpha', contentType=InputData.FloatType, strictMode=True))
    param.addSub(InputData.parameterInputFactory('A'    , contentType=InputData.FloatType, strictMode=True))

    inputSpecification.addSub(param)
    return inputSpecification

  def __init__(self):
    """
      Default Constructor
    """
    GradientBasedOptimizer.__init__(self)
    self.paramDict['pertSingleGrad'] = 1
    self.currentDirection = None
    self.stochasticDistribution = None                        # Distribution used to generate perturbations
    self.stochasticEngine = None                              # Random number generator used to generate perturbations
    self.stochasticEngineForConstraintHandling = Distributions.returnInstance('Normal',self)
    self.stochasticEngineForConstraintHandling.mean, self.stochasticEngineForConstraintHandling.sigma = 0, 1
    self.stochasticEngineForConstraintHandling.upperBoundUsed, self.stochasticEngineForConstraintHandling.lowerBoundUsed = False, False
    self.stochasticEngineForConstraintHandling.initializeDistribution()

  def localInputAndChecks(self, xmlNode):
    """
      Local method for additional reading.
      @ In, xmlNode, xml.etree.ElementTree.Element, Xml element node
      @ Out, None
    """
    GradientBasedOptimizer.localInputAndChecks(self, xmlNode)
    self.currentDirection   = None
    self.paramDict['alpha'] = float(self.paramDict.get('alpha', 0.602))
    self.paramDict['gamma'] = float(self.paramDict.get('gamma', 0.101))
    self.paramDict['A']     = float(self.paramDict.get('A', self.limit['mdlEval']/10.))
    self.paramDict['a']     = self.paramDict.get('a', None)
    self.paramDict['c']     = float(self.paramDict.get('c', 0.005))
    #FIXME the optimization parameters should probably all operate ONLY on normalized data!
    #  -> perhaps the whole optimizer should only work on optimized data.

    numValues = self._numberOfSamples()

    #FIXME normalizing doesn't seem to have the desired effect, currently; it makes the step size very small (for large scales)
    #if "a" was defaulted, use the average scale of the input space.
    #This is the suggested value from the paper, missing a 1/gradient term since we don't know it yet.
    if self.paramDict['a'] is None:
      self.paramDict['a'] = mathUtils.hyperdiagonal(np.ones(numValues)) # the features are always normalized
      self.raiseAMessage('Defaulting "a" gradient parameter to',self.paramDict['a'])
    else:
      self.paramDict['a'] = float(self.paramDict['a'])

    self.constraintHandlingPara['innerBisectionThreshold'] = float(self.paramDict.get('innerBisectionThreshold', 1e-2))
    self.constraintHandlingPara['innerLoopLimit'] = float(self.paramDict.get('innerLoopLimit', 1000))

    self.gradDict['pertNeeded'] = self.gradDict['numIterForAve'] * (self.paramDict['pertSingleGrad']+1)

    # determine the number of indpendent variables (scalar and vectors included)
    stochDist = self.paramDict.get('stochasticDistribution', 'Hypersphere')
    if stochDist == 'Bernoulli':
      self.stochasticDistribution = Distributions.returnInstance('Bernoulli',self)
      self.stochasticDistribution.p = 0.5
      self.stochasticDistribution.initializeDistribution()
      # Initialize bernoulli distribution for random perturbation. Add artificial noise to avoid that specular loss functions get false positive convergence
      # FIXME there has to be a better way to get two random numbers
      self.stochasticEngine = lambda: [(0.5+randomUtils.random()*(1.+randomUtils.random()/1000.*randomUtils.randomIntegers(-1, 1, self))) if self.stochasticDistribution.rvs() == 1 else
                                   -1.*(0.5+randomUtils.random()*(1.+randomUtils.random()/1000.*randomUtils.randomIntegers(-1, 1, self))) for _ in range(numValues)]
    elif stochDist == 'Hypersphere':
      self.stochasticEngine = lambda: randomUtils.randPointsOnHypersphere(numValues) if numValues > 1 else [randomUtils.randPointsOnHypersphere(numValues)]
    else:
      self.raiseAnError(IOError, self.paramDict['stochasticEngine']+'is currently not supported for SPSA')

  def localLocalInitialize(self, solutionExport):
    """
      Method to initialize local settings.
      @ In, solutionExport, DataObject, a PointSet to hold the solution
      @ Out, None
    """
    self._endJobRunnable = (self._endJobRunnable*self.gradDict['pertNeeded'])+len(self.optTraj)

  def _newOptPointAdd(self, gradient, traj):
    """
      This local method add a new opt point based on the gradient
      @ In, gradient, dict, dictionary containing the gradient
      @ In, traj, int, trajectory
      @ Out, None
    """
    ak = self._computeGainSequenceAk(self.paramDict,self.counter['varsUpdate'][traj],traj) # Compute the new ak
    self.optVarsHist[traj][self.counter['varsUpdate'][traj]] = {}
    varK = dict((var,self.counter['recentOptHist'][traj][0][var]) for var in self.getOptVars(traj))
    varKPlus,modded = self._generateVarsUpdateConstrained(traj,ak,gradient,varK)
    #check for redundant paths
    if len(self.optTrajLive) > 1 and self.counter['solutionUpdate'][traj] > 0:
      self._removeRedundantTraj(traj, varKPlus)
    #if the new point was modified by the constraint, reset the step size
    if modded:
      del self.counter['lastStepSize'][traj]
      self.raiseADebug('Resetting step size for trajectory',traj,'due to hitting constraints')
    self.queueUpOptPointRuns(traj,varKPlus)

  def localStillReady(self, ready, convergence = False):
    """
      Determines if optimizer is ready to provide another input.  If not, and if jobHandler is finished, this will end sampling.
      @ In, ready, bool, variable indicating whether the caller is prepared for another input.
      @ In, convergence, bool, optional, variable indicating whether the convergence criteria has been met.
      @ Out, ready, bool, variable indicating whether the caller is prepared for another input.
    """
    self.nextActionNeeded = (None,None) #prevents carrying over from previous run
    #get readiness from parent
    ready = ready and GradientBasedOptimizer.localStillReady(self,ready,convergence)
    #if not ready, just return that
    if not ready:
      return ready
    # loop over live trajectories and look for any actions that are ready to proceed
    # -> indicate the first ready action through the self.nextActionNeeded mechanic
    self.raiseADebug('Reviewing status of trajectories:')
    for traj in self.optTraj:
      self.raiseADebug('   Traj: "{:^n}": Process: "{:^30.30}", Reason: "{:^30.30}"'.format(traj,str(self.status[traj].get('process','None')),self.status[traj].get('reason','None'),n=len(str(max(self.optTraj)))))
    for _ in range(len(self.optTrajLive)):
      # despite several attempts, this is the most elegant solution I've found to assure each
      #   trajectory gets even treatment.
      traj = self.optTrajLive.pop(0)
      self.optTrajLive.append(traj)
      self.raiseADebug('Checking readiness for traj "{}" ...'.format(traj))
      process = self.status[traj]['process']
      reason = self.status[traj]['reason']

      # still in the process of submitting new points for evaluating the gradient
      if process == 'submitting grad eval points':
        self.nextActionNeeded = ('add new grad evaluation point',traj)
        break

      # still in the process of submitting new opt point evaluations (for stochastic denoising)
      elif process == 'submitting new opt points':
        if reason == 'just started':
          self.nextActionNeeded = ('start new trajectory',traj)
          break
        elif reason in ['seeking new opt point','received recommended point']:
          self.nextActionNeeded = ('add more opt point evaluations',traj)
          break
        elif reason in 'failed run':
          self.nextActionNeeded = ('add more opt point evaluations',traj)
          if len(self.submissionQueue[traj]) == 0:
            gradient = self.counter['gradientHistory'][traj][0]
            self._newOptPointAdd(gradient,traj)
        else:
          self.raiseAnError(RuntimeError,'unexpected reason for submitting new opt points:',reason)

      # all gradient evaluation points submitted, but waiting for them to be collected
      elif process == 'collecting grad eval points':
        # check to see if the grad evaluation points have all been collected
        evalsFinished = True
        for pertID in range(self.gradDict['pertNeeded']):
          if not self._checkModelFinish(traj,self.counter['varsUpdate'][traj],pertID)[0]:
            evalsFinished = False
            break
        # if grad eval pts are finished, then evaluate the gradient
        if evalsFinished:
          # collect output values for perturbed points
          #for i in range(1,self.gradDict['numIterForAve']*2,2):
          for i in self.perturbationIndices:
            evalIndex = self._checkModelFinish(traj,self.counter['varsUpdate'][traj],i)[1]
            outval = float(self.mdlEvalHist.realization(index=evalIndex)[self.objVar])
            self.gradDict['pertPoints'][traj][i]['output'] = outval
          # reset per-opt-point counters, forward the varsUpdate
          self.counter['perturbation'][traj] = 0
          self.counter['varsUpdate'][traj] += 1
          # evaluate the gradient
          gradient = self.evaluateGradient(self.gradDict['pertPoints'][traj],traj)
          # establish a new point, if found; FIXME otherwise?
          if len(self.submissionQueue[traj]) == 0:
            self._newOptPointAdd(gradient,traj)
            # if trajectory was killed for redundancy, continue on to check next trajectory for readiness
            if self.status[traj]['reason'] == 'removed as redundant':
              continue # loops back to the next opt traj
          self.nextActionNeeded = ('add more opt point evaluations',traj)
          self.status[traj]['process'] = 'submitting new opt points'
          self.status[traj]['reason'] = 'seeking new opt point'
          break
        # otherwise, we're not ready to sample yet
        else:
          self.raiseADebug('    Traj "{}": Waiting on collection of gradient evaluation points.'.format(traj))
          continue

      # all opt point evaluations have been submitted, but waiting for them to be collected
      elif process == 'collecting new opt points':
        self.raiseADebug('    Traj "{}": Waiting on collection of new optimization point.'.format(traj))
        continue

      # trajectory has already converged, so there's no work to do for it.
      elif reason == 'converged':
        self.raiseADebug('    Traj "{}": Trajectory is marked as converged.'.format(traj))
        continue

      # trajectory didn't converge, but was killed because of redundancy
      elif reason == 'removed as redundant':
        self.raiseADebug('    Traj "{}": Trajectory is removed for redundancy.'.format(traj))
        continue
      # unknown status
      else:
        self.raiseAnError(RuntimeError,'Unrecognized status:'.format(traj),self.status[traj])
    # end loop through trajectories looking for new actions

    # if we did not find an action, we're not ready to provide an input
    if self.nextActionNeeded[0] is None:
      self.raiseADebug('Not ready to provide a sample yet.')
      return False
    else:
      self.raiseADebug('Next action needed: "%s" on trajectory "%i"' %self.nextActionNeeded)
      return True

  def localEvaluateGradient(self, optVarsValues, traj,  gradient = None):
    """
      Local method to evaluate gradient.
      @ In, optVarsValues, dict, dictionary containing perturbed points.
                                 optVarsValues should have the form {pertIndex: {varName: [varValue1 varValue2]}}
                                 Therefore, each optVarsValues[pertIndex] should return a dict of variable values
                                 that is sufficient for gradient evaluation for at least one variable
                                 (depending on specific optimization algorithm)
      @ In, traj, int, the trajectory id
      @ In, gradient, dict, optional, dictionary containing gradient estimation by the caller.
                                      gradient should have the form {varName: gradEstimation}
      @ Out, gradient, dict, dictionary containing gradient estimation. gradient should have the form {varName: gradEstimation}
    """
    gradArray = {}
    # number of gradient evaluations to consider (denoising or resampling)
    numGrads = self.gradDict['numIterForAve']
    # prepopulate array of collected gradient
    for var in self.getOptVars(traj=traj):
      gradArray[var] = np.zeros(numGrads,dtype=object)
    # Evaluate gradient at each point
    for i in range(numGrads):
      opt  = optVarsValues[i]            #the latest opt point
      pert = optVarsValues[i + numGrads] #the perturbed point
      # calculate grad(F) wrt each input variable
      # fix infinities!
      lossDiff = mathUtils.diffWithInfinites(pert['output'],opt['output'])
      #cover "max" problems
      # TODO it would be good to cover this in the base class somehow, but in the previous implementation this
      #   sign flipping was only called when evaluating the gradient.
      #   Perhaps the sign should flip when evaluating the next point to take, instead of forcing gradient descent
      if self.optType == 'max':
        lossDiff *= -1.0
      for var in self.getOptVars(traj=traj):
        # NOTE: gradient is calculated in normalized space
        dh = pert['inputs'][var] - opt['inputs'][var]
        # a sample so close cannot be taken without violating minimum step, so this check should not be necessary (left for reference)
        #if abs(dh) < 1e-15:
        #  self.raiseAnError(RuntimeError,'While calculating the gradArray a "dh" of zero was found for var:',var)
        gradArray[var][i] = lossDiff/dh
    gradient = {}
    for var in self.getOptVars(traj=traj):
      mean = gradArray[var].mean()
      gradient[var] = np.atleast_1d(mean)
    return gradient

  def localGenerateInput(self,model,oldInput):
    """
      Method to generate input for model to run
      @ In, model, model instance, it is the instance of a RAVEN model
      @ In, oldInput, list, a list of the original needed inputs for the model (e.g. list of files, etc. etc)
      @ Out, None
    """
    GradientBasedOptimizer.localGenerateInput(self,model,oldInput)
    action, traj = self.nextActionNeeded
    #store traj as active for sampling
    self.inputInfo['trajID'] = traj+1
    self.inputInfo['varsUpdate'] = self.counter['varsUpdate'][traj]
    #"action" and "traj" are set in localStillReady
    #"action" is a string of the next action needed by the optimizer in order to move forward
    #"traj" is the trajectory that is in need of the action

    if action == 'start new trajectory':
      if len(self.submissionQueue[traj]) == 0:
        self.optVarsHist[traj][self.counter['varsUpdate'][traj]] = {}
        values = {}
        for var in self.getOptVars(traj=traj):
          values[var] = self.optVarsInit['initial'][var][traj]
          #if exceeding bounds, bring value within 1% of range
          values[var] = self._checkBoundariesAndModify(self.optVarsInit['upperBound'][var],
                                                            self.optVarsInit['lowerBound'][var],
                                                            self.optVarsInit['ranges'][var],values[var],0.99,0.01)
        # record this suggested new optimal point (in normalized input space)
        data = self.normalizeData(values)
        self.updateVariableHistory(values,traj)
        # set up a new batch of runs on the new optimal point (batch size 1 unless more requested by user)
        self.queueUpOptPointRuns(traj,data)
      #"submit" the next queued point
      prefix,point = self.getQueuedPoint(traj)
      for var in self.getOptVars(traj=traj):
        self.values[var] = point[var]
      self.inputInfo['prefix'] = prefix
      # if all iterations have been submitted, the optimizer is now in collection mode.
      if len(self.submissionQueue[traj]) == 0:
        self.status[traj]['process'] = 'collecting new opt points'

    elif action == 'add new grad evaluation point':
      #increment submission number
      self.counter['perturbation'][traj] += 1
      #if this is the first perturbation, prep all the perturbation (aka gradient evaluation) points we need to run
      #note that currently we use the opt point as half of the perturbation point, so each gradient eval will be between
      #   the opt point and a different perturbed point
      if self.counter['perturbation'][traj] == 1:
        # Generate all the perturbations at once, then we can submit them one at a time
        ck = self._computeGainSequenceCk(self.paramDict,self.counter['varsUpdate'][traj]+1)
        varK = dict((var,self.counter['recentOptHist'][traj][0][var]) for var in self.getOptVars(traj))
        #check the submission queue is empty; otherwise something went wrong # TODO this is a sanity check, might be removed for efficiency
        assert(len(self.submissionQueue[traj])==0)
        for i in self.perturbationIndices:
          direction = self._getPerturbationDirection(i, traj)
          point = {}
          index = 0
          for var in self.getOptVars(traj=traj):
            size = np.prod(self.variableShapes[var])
            if size > 1:
              new = np.zeros(size)
              for v,origVal in enumerate(varK[var]):
                new[v] = origVal + ck*direction[index]
                new[v] = self._checkBoundariesAndModify(1.0, 0.0, 1.0, new[v], 0.9999, 0.0001)
                index += 1
              point[var] = new
            else:
              val = varK[var] + ck*direction[index]
              index += 1
              val = self._checkBoundariesAndModify(1.0, 0.0, 1.0, val, 0.9999, 0.0001)
              point[var] = val
          #create identifier
          prefix = self._createEvaluationIdentifier(traj,self.counter['varsUpdate'][traj],i)
          #queue it up
          self.submissionQueue[traj].append({'inputs':point,'prefix':prefix})
      #end if-first-time conditional
      #get a queued entry to run
      entry = self.submissionQueue[traj].popleft()
      prefix = entry['prefix']
      point = entry['inputs']
      self.gradDict['pertPoints'][traj][int(prefix.split('_')[-1])] = {'inputs':point}#self.normalizeData(point)}
      point = self.denormalizeData(point)
      for var in self.getOptVars(traj=traj):
        self.values[var] = point[var]
      # use 'prefix' to locate the input sent out, as <traj>_<varUpdate>_<eval number>, as 0_0_2
      self.inputInfo['prefix'] = prefix
      # if all required points are submitted, switch into collection mode
      if len(self.submissionQueue[traj]) == 0:
        self.status[traj]['process'] = 'collecting grad eval points'

    elif action == 'add more opt point evaluations':
      #take a sample from the queue
      prefix,point = self.getQueuedPoint(traj)
      #prep the point for running
      for var in self.getOptVars(traj=traj):
        self.values[var] = point[var]
      self.updateVariableHistory(self.values,traj)
      # use 'prefix' to locate the input sent out
      self.inputInfo['prefix'] = prefix
      # if all points submitted (after current one gets submitted), switch to collection mode
      if len(self.submissionQueue[traj]) == 0:
        self.status[traj]['process'] = 'collecting new opt points'

    #unrecognized action
    else:
      self.raiseAnError(RuntimeError,'Unrecognized "action" in localGenerateInput:',action)
    self.raiseADebug('Queuing run "{}"'.format(self.inputInfo['prefix']))

  def _checkBoundariesAndModify(self,upperBound,lowerBound,varRange,currentValue,pertUp,pertLow):
    """
      Method to check the boundaries and add a perturbation in case they are violated
      @ In, upperBound, float, the upper bound for the variable
      @ In, lowerBound, float, the lower bound for the variable
      @ In, varRange, float, the variable range
      @ In, currentValue, float, the current value
      @ In, pertUp, float, the perturbation to apply in case the upper bound is violated
      @ In, pertLow, float, the perturbation to apply in case the lower bound is violated
      @ Out, convertedValue, float, the modified value in case the boundaries are violated
    """
    # if variable is a vector of values, treat each independently
    if isinstance(currentValue,np.ndarray):
      convertedValue = np.asarray([self._checkBoundariesAndModify(upperBound,lowerBound,varRange,value,pertUp,pertLow) for value in currentValue])
    else:
      convertedValue = currentValue
      if currentValue > upperBound:
        convertedValue = pertUp*varRange + lowerBound
      elif currentValue < lowerBound:
        convertedValue = pertLow*varRange + lowerBound
    return convertedValue

  def clearCurrentOptimizationEffort(self,traj):
    """
      See base class.  Used to clear out current optimization information and start clean.
      For the SPSA, this means clearing out the perturbation points
      @ In, traj, int, index of trajectory being cleared
      @ Out, None
    """
    self.raiseADebug('Clearing current optimization efforts ...')
    self.counter ['perturbation'   ][traj] = 0
    self.counter ['gradientHistory'][traj] = [{},{}]
    self.counter ['gradNormHistory'][traj] = [0,0]
    #only clear non-opt points from pertPoints
    for i in self.perturbationIndices:
      self.gradDict['pertPoints'][traj][i] = 0
    self.convergeTraj[traj] = False
    self.status[traj] = {'process':'submitting grad eval points','reason':'found new opt point'}
    try:
      del self.counter['lastStepSize'][traj]
    except KeyError:
      pass
    try:
      del self.recommendToGain[traj]
    except KeyError:
      pass

  def _updateParameters(self):
    """
      Uses information about the gradient and optimizer trajectory to update parameters
      Updated parameters include [a,A,alpha,c,gamma]
    """
    pass #future tool
    #update A <-- shouldn't be needed since A is an early-life stabilizing parameter
    #update alpha <-- if we're moving basically in the same direction, don't damp!
    #update a <-- increase or decrease step size based on gradient information
    #  determine the minimum desirable step size at early calcualtion (1/10 of average
    #update c <-- estimate stochasticity of the response; if low, "c" tends to 0
    #update gamma <-- distance between samples to determine gradient.  Scales with step size?

  def _generateVarsUpdateConstrained(self,traj,ak,gradient,varK):
    """
      Method to generate input for model to run, considering also that the input satisfies the constraint
      @ In, traj, int, trajectory label for whom we are generating variables with constraint consideration
      @ In, ak, float or array, it is gain for variable update (if array, different gain for each variable)
      @ In, gradient, dictionary, contains the gradient information for variable update
      @ In, varK, dictionary, current variable values (normalized)
      @ Out, varKPlus, dictionary, variable values for next iteration.
      @ Out, modded, bool, if True the point was modified by the constraint
    """
    varKPlus = {}
    gain = [ak]*self._numberOfSamples() #technically too many entries, but unneeded ones will be *0 anyway just below here
    gain = np.asarray(gain)
    index = 0
    for var in self.getOptVars(traj=traj):
      numSamples = np.prod(self.variableShapes[var])
      if numSamples == 1:
        new = varK[var]-gain[index]*gradient.get(var,0.0)
        index += 1
      else:
        new = np.zeros(numSamples)
        for i in range(numSamples):
          new[i] = varK[var][i] - gain[index] * gradient.get(var,[0.0]*i)[i]
          index += 1
      varKPlus[var] = new
    satisfied, activeViolations = self.checkConstraint(self.denormalizeData(varKPlus))
    if satisfied:
      return varKPlus, False
    # else if not satisfied ...
    # check if the active constraints are the boundary ones. In this case, try to project the gradient at an angle
    modded = False
    if len(activeViolations['internal']) > 0:
      self.raiseADebug('Attempting to fix constraint violation with gradient projection ...')
      modded = True
      projectedOnBoundary= {}
      for var,under,over in activeViolations['internal']:
        if np.prod(self.variableShapes[var]) == 1:
          if np.sum(over) > 0:
            projectedOnBoundary[var] = self.optVarsInit['upperBound'][var]
          elif np.sum(under) > 0:
            projectedOnBoundary[var] = self.optVarsInit['lowerBound'][var]
          gradient[var] = 0.0
        else:
          projectedOnBoundary[var] = self.denormalizeData({var:varKPlus[var]})[var]
          projectedOnBoundary[var][under] = self.optVarsInit['lowerBound'][var]
          projectedOnBoundary[var][over] = self.optVarsInit['upperBound'][var]
          gradient[var][np.logical_or(under,over)] = 0.0
      varKPlus.update(self.normalizeData(projectedOnBoundary))
      newNormWithoutComponents = self.calculateMultivectorMagnitude(gradient.values())
      for var in gradient.keys():
        gradient[var] = gradient[var]/newNormWithoutComponents if newNormWithoutComponents != 0.0 else gradient[var]

    if len(activeViolations['external']) == 0:
      return varKPlus, modded

    self.raiseADebug('Attempting to fix constraint violation by shortening gradient vector ...')
    # Try to find varKPlus by shorten the gradient vector
    self.raiseADebug('Trajectory "{}" hit constraints ...'.format(traj))
    self.raiseADebug('  Attempting to shorten step length ...')
    foundVarsUpdate, varKPlus = self._bisectionForConstrainedInput(traj,varK, ak, gradient)
    if foundVarsUpdate:
      self.raiseADebug('   ... successfully found new point by shortening length.')
      return varKPlus, True

    self.raiseADebug('Attempting to fix constraint violation by rotating towards orthogonal ...')
    # Try to find varKPlus by rotate the gradient towards its orthogonal, since we consider the gradient as perpendicular
    # with respect to the constraints hyper-surface
    self.raiseADebug('  Attempting instead to rotate trajectory ...')
    innerLoopLimit = self.constraintHandlingPara['innerLoopLimit']
    if innerLoopLimit < 0:
      self.raiseAnError(IOError, 'Limit for internal loop for constraint handling should be nonnegative')
    loopCounter = 0
    foundPendVector = False
    # search for the perpendicular vector
    while not foundPendVector and loopCounter < innerLoopLimit:
      loopCounter += 1
      # randomly choose the index of a variable to be the dependent? pivot
      depVarPos = randomUtils.randomIntegers(0,len(self.getOptVars(traj=traj))-1,self)
      # if that variable is multidimensional, pick a dimension -> this is not precisely equal probability of picking, but that should be okay.
      varSize = np.prod(self.variableShapes[var])
      if varSize > 1:
        depVarIdx = randomUtils.randomIntegers(0,varSize-1,self)
      pendVector = {}
      npDot = 0
      for varID, var in enumerate(self.getOptVars(traj=traj)):
        varSize = np.prod(self.variableShapes[var])
        if varSize == 1:
          pendVector[var] = self.stochasticEngineForConstraintHandling.rvs() if varID != depVarPos else 0.0
          npDot += pendVector[var]*gradient[var]
        else:
          for i in range(varSize):
            pendVector[var][i] = self.stochasticEngineForConstraintHandling.rvs() if (varID != depVarPos and depVarIdx != i) else 0.0
          npDot += np.sum(pendVector[var]*gradient[var])
      # TODO does this need to be in a separate loop or can it go with above?
      for varID, var in enumerate(self.getOptVars(traj=traj)):
        if varID == depVarPos:
          varSize = np.prod(self.variableShapes[var])
          if varSize == 1:
            pendVector[var] = -npDot/gradient[var]
          else:
            pendVector[var][depVarIdx] = -npDot/gradient[var][depVarIdx]

      r  = self.calculateMultivectorMagnitude([  gradient[var] for var in self.getOptVars(traj=traj)])
      r /= self.calculateMultivectorMagnitude([pendVector[var] for var in self.getOptVars(traj=traj)])
      for var in self.getOptVars(traj=traj):
        pendVector[var] = copy.deepcopy(pendVector[var])*r

      varKPlus = {}
      index = 0
      for var in self.getOptVars(traj=traj):
        varSize = np.prod(self.variableShapes[var])
        new = np.zeros(varSize)
        for i in range(varSize):
          if varSize == 1:
            new = copy.copy(varK[var]-gain[index]*pendVector[var]*1.0)
          else:
            new[i] = copy.copy(varK[var][i]-gain[index]*pendVector[var][i]*1.0)
          index += 1
        varKPlus[var] = new

      foundPendVector, activeConstraints = self.checkConstraint(self.denormalizeData(varKPlus))
      if not foundPendVector:
        foundPendVector, varKPlus = self._bisectionForConstrainedInput(traj,varK, gain, pendVector)
      gain = gain/2.

    if foundPendVector:
      lenPendVector = 0
      for var in self.getOptVars(traj=traj):
        lenPendVector += np.sum(pendVector[var]**2)
      lenPendVector = np.sqrt(lenPendVector)

      rotateDegreeUpperLimit = 2
      while self.angleBetween(traj,gradient, pendVector) > rotateDegreeUpperLimit:
        sumVector = {}
        lenSumVector = 0
        for var in self.getOptVars(traj=traj):
          sumVector[var] = gradient[var] + pendVector[var]
          lenSumVector += np.sum(sumVector[var]**2)

        tempTempVarKPlus = {}
        index = 0
        for var in self.getOptVars(traj=traj):
          sumVector[var] = copy.deepcopy(sumVector[var]/np.sqrt(lenSumVector)*lenPendVector)
          varSize = np.prod(self.variableShapes[var])
          new = np.zeros(varSize)
          for i in range(varSize):
            if varSize == 1:
              new = copy.copy(varK[var]-gain[index]*sumVector[var]*1.0)
            else:
              new[i] = copy.copy(varK[var][i]-gain[index]*sumVector[var][i]*1.0)
            index += 1
          tempTempVarKPlus[var] = new
        satisfied, activeConstraints = self.checkConstraint(self.denormalizeData(tempTempVarKPlus))
        if satisfied:
          varKPlus = copy.deepcopy(tempTempVarKPlus)
          pendVector = copy.deepcopy(sumVector)
        else:
          gradient = copy.deepcopy(sumVector)
      self.raiseADebug('   ... successfully found new point by rotating trajectory.')
      return varKPlus, True
    varKPlus = varK
    self.raiseADebug('   ... did not successfully find new point.')
    return varKPlus, False

  def _bisectionForConstrainedInput(self,traj,varK,ak,vector):
    """
      Method to find the maximum fraction of the step size "ak" such that taking a step of size "ak" in
      the direction of "vector" starting from current input point "varK" will not violate the constraint function.
      @ In, traj, int, trajectory label for whom we are working
      @ In, varK, dictionary, current variable values
      @ In, ak, float or array, it is step size (gain) for variable update (if array, different gain for each variable)
      @ In, vector, dictionary, contains the gradient information as a unit vector
      @ Out, _bisectionForConstrainedInput, tuple(bool,dict), (indicating whether a fractional step is found, contains the fraction of step that satisfies constraint)
    """
    try:
      gain = ak[:]
    except (TypeError,IndexError):
      gain = [ak]*self._numberOfSamples() #technically incorrect, but missing ones will be *0 anyway just below here

    innerBisectionThreshold = self.constraintHandlingPara['innerBisectionThreshold']
    if innerBisectionThreshold <= 0 or innerBisectionThreshold >= 1: #FIXME REWORK this is an input check, not a runtime check
      self.raiseAnError(ValueError, 'The innerBisectionThreshold should be greater than 0 and less than 1')
    bounds = [0, 1.0]
    tempVarNew = {}
    frac = 0.5
    while np.absolute(bounds[1]-bounds[0]) >= innerBisectionThreshold:
      index = 0
      for var in self.getOptVars(traj=traj):
        numSamples = np.prod(self.variableShapes[var])
        new = np.zeros(numSamples)
        for i in range(numSamples):
          if numSamples > 1:
            new[i] = copy.copy(varK[var][i] - gain[index] * vector[var][i]*1.0*frac) # FIXME is this copy needed?
          else:
            new = copy.copy(varK[var]-gain[index]*vector[var]*1.0*frac) # FIXME is this copy needed?
          index += 1
        tempVarNew[var] = new

      satisfied,_ = self.checkConstraint(tempVarNew)
      if satisfied:
        bounds[0] = copy.deepcopy(frac)
        if np.absolute(bounds[1]-bounds[0]) >= innerBisectionThreshold:
          varKPlus = copy.deepcopy(tempVarNew)
          return True, varKPlus
        frac = copy.deepcopy(bounds[1]+bounds[0])/2.0
      else:
        bounds[1] = copy.deepcopy(frac)
        frac = copy.deepcopy(bounds[1]+bounds[0])/2.0
    return False, None

  def angleBetween(self, traj, d1, d2):
    """
      Evaluate the angle between the two dictionaries of vars (d1 and d2) by means of the dot product. Unit: degree
      @ In, traj, int, trajectory label for whom we are working
      @ In, d1, dict, first vector
      @ In, d2, dict, second vector
      @ Out, angleD, float, angle between d1 and d2 with unit of degree
    """
    nVar = len(self.getOptVars())
    v1, v2 = np.zeros(shape=[nVar,]), np.zeros(shape=[nVar,])
    for cnt, var in enumerate(self.getOptVars(traj=traj)):
      v1[cnt], v2[cnt] = copy.deepcopy(d1[var]), copy.deepcopy(d2[var])
    angle = np.arccos(np.dot(v1, v2)/np.linalg.norm(v1)/np.linalg.norm(v2))
    if np.isnan(angle):
      if (v1 == v2).all():
        angle = 0.0
      else:
        angle = np.pi
    angleD = np.rad2deg(angle)
    return angleD

  def _computeGainSequenceCk(self,paramDict,iterNum):
    """
      Utility function to compute the ck coefficients (gain sequence ck)
      @ In, paramDict, dict, dictionary containing information to compute gain parameter
      @ In, iterNum, int, current iteration index
      @ Out, ck, float, current value for gain ck
    """
    c, gamma = paramDict['c'], paramDict['gamma']
    ck = c / (iterNum) ** gamma *1.0
    return ck

  def _computeGainSequenceAk(self,paramDict,iterNum,traj):
    """
      Utility function to compute the ak coefficients (gain sequence ak)
      @ In, paramDict, dict, dictionary containing information to compute gain parameter
      @ In, iterNum, int, current iteration index
      @ Out, new, float, current value for gain ak
    """
    #TODO FIXME is this a good idea?
    try:
      ak = self.counter['lastStepSize'][traj]
    except KeyError:
      a, A, alpha = paramDict['a'], paramDict['A'], paramDict['alpha']
      ak = a / (iterNum + A) ** alpha
    # modify step size based on the history of the gradients used
    frac = self.fractionalStepChangeFromGradHistory(traj)
    new = ak*frac
    self.raiseADebug('step gain size for traj "{}" iternum "{}": {:1.3e} (root {:1.2e} frac {:1.2e})'.format(traj,iterNum,new,ak,frac))
    self.counter['lastStepSize'][traj] = new
    return new

  def _getAlgorithmState(self,traj):
    """
      Returns values specific to this algorithm such that it could pick up again relatively easily from here.
      #REWORK functionally a getstate/setstate for multilevel
      @ In, traj, int, the trajectory being saved
      @ Out, state, dict, keys:values this algorithm cares about saving for this trajectory
    """
    state = {}
    state['lastStepSize']    = copy.deepcopy(self.counter['lastStepSize'   ].get(traj,None))
    state['gradientHistory'] = copy.deepcopy(self.counter['gradientHistory'].get(traj,None))
    state['recommendToGain'] = copy.deepcopy(self.recommendToGain           .get(traj,None))
    return state

  def _setAlgorithmState(self,traj,state):
    """
      #REWORK functionally a getstate/setstate for multilevel
      @ In, traj, int, the trajectory being saved
      @ In, state, dict, keys:values this algorithm cares about saving for this trajectory
      @ Out, None
    """
    if state is None:
      return
    if state['lastStepSize'] is not None:
      self.counter['lastStepSize'][traj] = state['lastStepSize']
    if state['gradientHistory'] is not None:
      self.counter['gradientHistory'][traj] = state['gradientHistory']
    if state['recommendToGain'] is not None:
      self.recommendToGain[traj] = state['recommendToGain']

  def _getPerturbationDirection(self,perturbationIndex, traj):
    """
      This method is aimed to get the perturbation direction (i.e. in this case the random perturbation versor)
      @ In, perturbationIndex, int, the perturbation index (stored in self.perturbationIndices)
      @ In, traj, int, the trajectory id
      @ Out, direction, list, the versor for each optimization dimension
    """
    if perturbationIndex == self.perturbationIndices[0]:
      direction = self.stochasticEngine()
      self.currentDirection = direction
    else:
      # in order to perform the de-noising we keep the same perturbation direction and we repeat the evaluation multiple times
      direction = self.currentDirection
    return direction

max_stars_repo_name: bonifak/raven

max_issues_repo_name: bonifak/raven

max_forks_repo_name: bonifak/raven


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

plt.style.use("zzo_color")
def get_data(path, x, y):
    with open(path) as f:
        lines = f.readlines()
    for line in lines:
        t, cwnd = line.split("\t")
        x.append(float(t))
        y.append(float(cwnd.strip()))

def plot_data(path, scale=1):
    x = []
    y = []
    get_data(path, x, y)
    plt.plot(x, np.array(y) * scale)

def main(k):
    x = []
    y = []

    path_cwnd = "/tmp/source{}_cwnd_data.tsv"
    path_rtt = "/tmp/source{}_rtt_data.tsv"

    legend = []
    for i in range(1, k+1):
        plot_data(path_cwnd.format(i))
        legend.append("$w_{{{}}}(t)$")
    bdp = 420
    plt.plot([0, 30], [420 / k, 420 / k], color="black", linestyle="dotted")
    legend.append("BDP")
    plt.legend(legend)
    plt.xlabel("time $t$ in \\si{\\second}")
    plt.ylabel("congestion window $w(t)$ in segments")
    tikzplotlib.save("/tmp/cwnd.tikz",
                     figurewidth=r"\textwidth", figureheight=r".25\textheight")

    plt.figure()
    for i in range(1, k+1):
        plot_data(path_rtt.format(i), 1e3)
    plt.xlabel("time $t$ in \\si{\\second}")
    plt.ylabel("RTT $\\tau(t)$ in \\si{\\milli\\second}")
    plt.yticks(range(140, 211, 10))
    tikzplotlib.save("/tmp/rtt.tikz",
                     figurewidth=r"\textwidth", figureheight=r".25\textheight")
    plt.show()

if __name__ == '__main__':
    try:
        k = int(sys.argv[1])
        assert k > 0
    except (IndexError, ValueError, AssertionError):
        print("Usage: python plot_graphs k")
        print("k is the number of flows that were emulated, k > 0!")
        exit(-1)
    main(k)

ext: py

max_stars_repo_name: GianGian/SDN-Slicing-in-ComNetsEmu

max_issues_repo_name: GianGian/SDN-Slicing-in-ComNetsEmu

max_forks_repo_name: hakiri/comnetsemu

reason: import numpy


text:
# Copyright 2016 - Alcatel-Lucent
#
# Licensed under the Apache License, Version 2.0 (the "License"); you may
# not use this file except in compliance with the License. You may obtain
# a copy of the License at
#
#      http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS, WITHOUT
# WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the
# License for the specific language governing permissions and limitations
# under the License.

from networkx.algorithms import components
from networkx.algorithms import simple_paths

from oslo_log import log as logging

from vitrage.common.constants import EdgeProperties as EProps
from vitrage.graph.algo_driver.algorithm import GraphAlgorithm
from vitrage.graph.algo_driver.algorithm import Mapping
from vitrage.graph.algo_driver.sub_graph_matching import NEG_CONDITION
from vitrage.graph.algo_driver.sub_graph_matching import subgraph_matching
from vitrage.graph.driver import Direction
from vitrage.graph.driver import Edge
from vitrage.graph.driver import Vertex
from vitrage.graph.filter import check_filter
from vitrage.graph.query import create_predicate

LOG = logging.getLogger(__name__)

class NXAlgorithm(GraphAlgorithm):

    def __init__(self, graph):
        """Create a new GraphAlgorithm

        :param graph: graph instance
        :type graph: driver.Graph
        """
        super(NXAlgorithm, self).__init__(graph)

    @classmethod
    def _create_new_graph(cls, *args, **kwargs):
        from vitrage.graph.driver.networkx_graph import NXGraph
        return NXGraph(args, **kwargs)

    def graph_query_vertices(self,
                             query_dict=None,
                             root_id=None,
                             depth=None,
                             direction=Direction.BOTH,
                             edge_query_dict=None):
        graph = self._create_new_graph('graph')

        if not root_id:
            root_id = self.graph.root_id
        root_data = self.graph._g.node[root_id]

        match_func = create_predicate(query_dict) if query_dict else None
        edge_match_func = create_predicate(edge_query_dict) \
            if edge_query_dict else None

        if match_func and not match_func(root_data):
            LOG.info('graph_query_vertices: root %s does not match filter %s',
                     str(root_id), str(query_dict))
            return graph

        n_result = []
        visited_nodes = set()
        n_result.append((root_id, self.graph.get_vertex(root_id).properties))
        e_result = []
        nodes_q = [(root_id, 0)]
        while nodes_q:
            node_id, curr_depth = nodes_q.pop(0)
            if (node_id in visited_nodes) or (depth and curr_depth >= depth):
                continue
            visited_nodes.add(node_id)
            (n_list, e_list) = self.graph._neighboring_nodes_edges_query(
                node_id,
                direction=direction,
                vertex_predicate=match_func,
                edge_predicate=edge_match_func)
            n_result.extend(n_list)
            e_result.extend(e_list)
            nodes_q.extend([(v_id, curr_depth + 1) for v_id, data in n_list])

        graph = self._create_new_graph(
            graph.name,
            graph.root_id,
            vertices=self._vertex_result_to_list(n_result),
            edges=self._edge_result_to_list(e_result))

        LOG.debug('graph_query_vertices: find graph: nodes %s, edges %s',
                  str(graph._g.nodes(data=True)),
                  str(graph._g.edges(data=True)))
        LOG.debug('graph_query_vertices: real graph: nodes %s, edges %s',
                  str(self.graph._g.nodes(data=True)),
                  str(self.graph._g.edges(data=True)))
        return graph

    def sub_graph_matching(self,
                           subgraph,
                           known_match,
                           validate=False):
        """Finds all the matching subgraphs in the graph

        In case the known_match has a subgraph edge with property
        "negative_condition" then run subgraph matching on the edge vertices
        and unite the results.
        Otherwise just run subgraph matching and return its result.

        :param subgraph: the subgraph to match
        :param known_match: starting point at the subgraph and the graph
        :param validate:
        :return: all the matching subgraphs in the graph
        """
        sge = known_match.subgraph_element
        ge = known_match.graph_element

        if not known_match.is_vertex and sge.get(NEG_CONDITION):
            source_matches = self._filtered_subgraph_matching(ge.source_id,
                                                              sge.source_id,
                                                              subgraph,
                                                              validate)
            target_matches = self._filtered_subgraph_matching(ge.target_id,
                                                              sge.target_id,
                                                              subgraph,
                                                              validate)

            return self._list_union(source_matches, target_matches)
        else:
            return subgraph_matching(self.graph,
                                     subgraph,
                                     [known_match],
                                     validate)

    def create_graph_from_matching_vertices(self,
                                            vertex_attr_filter=None,
                                            query_dict=None,
                                            edge_attr_filter=None):
        if query_dict:
            vertices = self.graph.get_vertices(query_dict=query_dict)
        elif vertex_attr_filter:
            vertices = self.graph.get_vertices(
                vertex_attr_filter=vertex_attr_filter)
        else:
            vertices = self.graph.get_vertices()

        vertices_ids = [vertex.vertex_id for vertex in vertices]

        graph = self._create_new_graph('graph')
        graph._g = self.graph._g.subgraph(vertices_ids)

        # delete non matching edges
        if edge_attr_filter:
            self._apply_edge_attr_filter(graph, edge_attr_filter)

        LOG.debug('match query, find graph: nodes %s, edges %s',
                  str(graph._g.nodes(data=True)),
                  str(graph._g.edges(data=True)))
        LOG.debug('match query, real graph: nodes %s, edges %s',
                  str(self.graph._g.nodes(data=True)),
                  str(self.graph._g.edges(data=True)))

        return graph

    def subgraph(self, entities):
        subgraph = self._create_new_graph('graph')
        subgraph._g = self.graph._g.subgraph(entities)
        return subgraph

    def connected_component_subgraphs(self, subgraph):
        return components.connected_component_subgraphs(
            subgraph._g.to_undirected(), copy=False)

    def all_simple_paths(self, source, target):
        return simple_paths.all_simple_paths(self.graph._g,
                                             source=source,
                                             target=target)

    def _filtered_subgraph_matching(self,
                                    ge_v_id,
                                    sge_v_id,
                                    subgraph,
                                    validate):
        """Runs subgraph_matching on edges vertices with filtering

        Runs subgraph_matching on edges vertices after checking if that vertex
        has real neighbors in the entity graph.
        """
        if self.graph.neighbors(ge_v_id,
                                edge_attr_filter={EProps.VITRAGE_IS_DELETED:
                                                  False}):
            template_vertex = subgraph.get_vertex(sge_v_id)
            graph_vertex = self.graph.get_vertex(ge_v_id)
            match = Mapping(template_vertex, graph_vertex, True)
            return subgraph_matching(self.graph, subgraph, [match], validate)

        return []

    @staticmethod
    def _edge_result_to_list(edge_result):
        d = dict()
        for source_id, target_id, label, data in edge_result:
            d[(source_id, target_id, label)] = \
                Edge(source_id, target_id, label, properties=data)
        return d.values()

    @staticmethod
    def _vertex_result_to_list(vertex_result):
        d = dict()
        for v_id, data in vertex_result:
            d[v_id] = Vertex(vertex_id=v_id, properties=data)
        return d.values()

    @staticmethod
    def _list_union(list_1, list_2):
        """Union of list that aren't hashable

        Can't use here set union because the items in the lists are
        dictionaries and they are not hashable for set.

        :return: list - union list
        """

        for target_item in list_2:
            if target_item not in list_1:
                list_1.append(target_item)

        return list_1

    @staticmethod
    def _apply_edge_attr_filter(graph, edge_attr_filter):
        for source, target, edge_data in graph._g.edges_iter(data=True):
            if not check_filter(edge_data, edge_attr_filter):
                graph._g.remove_edge(u=source, v=target)

reason: from networkx


text:
"""
This file is part of Data Analysis Recipes.
Copyright 2011, 2012 David W. Hogg (NYU).

Code to make a toy model of morphological star-galaxy separation
"""

import numpy as np
# this rc block must be before the matplotlib import?
from matplotlib import rc
rc('font',**{'family':'serif','serif':'Computer Modern Roman','size':18})
rc('text', usetex=True)
# now import matplotlib
import matplotlib
matplotlib.use('Agg')
import pylab as plt
import scipy.optimize as op

prefix = 'toy-morph'
suffix = 'png'
size_sigma = 0.05

def get_integral(alpha, S):
    return S**(1. - alpha) / (1. - alpha)

def integrate_over_bin(alpha, Sa, Sb):
    return get_integral(alpha, Sa) - get_integral(alpha, Sb)

def make_flux(alpha, Shi, Slo, size=(1)):
    return ((1. - alpha) * (np.random.uniform(size=size) * integrate_over_bin(alpha, Shi, Slo) + get_integral(alpha, Slo)))**(1. / (1. - alpha))

def observed_size(truesize):
    return np.sqrt(truesize**2 + 1.0) + size_sigma * np.random.normal(size=truesize.shape)

def mag_from_flux(flux):
    return 25. - 2.5 * np.log(flux) / np.log(10.)

def make_truth():
    Slo = 1.0
    Shi = 100.0
    nstar = 100
    ngala = 1000
    alphastar = 1.1
    alphagala = 1.7
    fluxstar = make_flux(alphastar, Slo, Shi, size=nstar)
    fluxgala = make_flux(alphagala, Slo, Shi, size=ngala)
    magstar = mag_from_flux(fluxstar)
    maggala = mag_from_flux(fluxgala)
    truesizestar = np.zeros_like(magstar)
    truesizegala = 0.4 * (25.5 - maggala) * np.random.gamma(3., size=maggala.shape)
    sizestar = observed_size(truesizestar)
    sizegala = observed_size(truesizegala)
    return magstar, sizestar, maggala, sizegala

def quantiles(x):
    nx = len(x)
    I = np.argsort(x)
    q = np.zeros(nx).astype(int)
    for i in range(4):
        q[I[(i * nx / 4):((i + 1) * nx / 4)]] = i
    splits = np.zeros(3).astype(float)
    for i in range(3):
        splits[i] = 0.5 * (np.max(x[q==i]) + np.min(x[q==(i+1)]))
    return q, splits

def straight_cut(t, c, ntotal=None):
    nc = len(c)
    cuts = np.arange(0.9,1.4,0.01)
    if ntotal is None:
        ntotal = np.sum(c == 0)
    completeness = np.array([float(np.sum(c[t < cut] == 0)) / float(ntotal) for cut in cuts])
    purity = np.array([float(np.sum(c[t < cut] == 0)) / float(np.sum(c[t < cut] < 3)) for cut in cuts])
    return cuts, completeness, purity

def prob_cut(logr, c):
    nc = len(c)
    cuts = np.arange(-5., 5., 0.1)
    completeness = np.array([float(np.sum(c[logr > cut] == 0)) / float(np.sum(c == 0)) for cut in cuts])
    purity = np.array([float(np.sum(c[logr > cut] == 0)) / float(np.sum(c[logr > cut] < 3)) for cut in cuts])
    return cuts, completeness, purity

def Gaussian(x, m, sigma):
    return np.exp(-0.5 * (x - m)**2 / sigma**2) / np.sqrt(2. * np.pi * sigma**2)

def likelihood_star(t):
    return Gaussian(t, 1.0, size_sigma)

np.random.seed(21)
gamma_sample = np.random.gamma(3., size=20000)
gamma_weight = np.ones_like(gamma_sample)
gamma_weight /= np.sum(gamma_weight)

# assumes the t, m come in one gala at a time
# scale is the scale of the size distribution, not of an individual galaxy
def likelihood_gala(ts, scale):
    return np.array([np.sum(gamma_weight * Gaussian(t, np.sqrt(1. + (scale * gamma_sample)**2), size_sigma)) for t in ts])

def probability(ts, pstar, scale):
    lstar = likelihood_star(ts)
    lgala = likelihood_gala(ts, scale)
    return pstar * lstar + (1. - pstar) * lgala

def total_log_probability(ts, pstar, scale):
    return np.sum(np.log(probability(ts, pstar, scale)))

def cost(pars, args):
    logodds, scale = pars
    pstar = np.exp(logodds) / (1. + np.exp(logodds))
    print pstar, scale
    ts = args
    return -1. * total_log_probability(ts, pstar, scale)

def main():
    np.random.seed(42)

    mstar, tstar, mgala, tgala = make_truth()
    cstar = np.zeros_like(mstar).astype(int)
    cgala = np.ones_like(mgala).astype(int)
    m = np.append(mstar, mgala)
    t = np.append(tstar, tgala)
    c = np.append(cstar, cgala)

    plt.clf()
    plt.hist(mgala, bins=5, histtype='step', color='b')
    plt.hist(mstar, bins=5, histtype='step', color='g')
    plt.xlim(20., 25.)
    plt.xlabel('magnitude $m$')
    plt.semilogy()
    plt.savefig('%s-hist.%s' % (prefix, suffix))

    plt.clf()
    plt.plot(tgala, mgala, 'bo', mew=0, alpha=0.5)
    plt.plot(tstar, mstar, 'go', mew=0, alpha=0.5)
    sizelim = (0.0, 10.0)
    plt.xlim(sizelim)
    plt.xlabel(r'size $\theta$')
    maglim = (25., 20.)
    plt.ylim(maglim)
    plt.ylabel('magnitude $m$')
    plt.savefig('%s-labeled-data.%s' % (prefix, suffix))

    plt.clf()
    plt.plot(t, m, 'ko', mew=0, alpha=0.5)
    plt.xlim(sizelim)
    plt.xlabel(r'size $\theta$')
    plt.ylim(maglim)
    plt.ylabel('magnitude $m$')
    plt.savefig('%s-data.%s' % (prefix, suffix))
    quants, splits = quantiles(m)
    for split in splits:
        plt.axhline(split, color='r', lw=2., alpha=0.75)
    plt.savefig('%s-data-qs.%s' % (prefix, suffix))

    lstar = likelihood_star(t)
    lgala = np.zeros_like(lstar)
    pstar = np.zeros_like(lstar)
    pgala = np.zeros_like(lstar)

    pstarbest = np.array([])
    scalebest = np.array([])
    lo, scale = 0.0, 1.0
    for q in range(4):
        I = (quants == q)
        tfit = t[I]
        lo, scale = op.fmin(cost, (lo, scale), args=(tfit,))
        thispstar = np.exp(lo) / (1. + np.exp(lo))
        pstarbest = np.append(pstarbest, thispstar)
        scalebest = np.append(scalebest, scale)
        lgala[I] = likelihood_gala(t[I], scale)
        pstar[I] = thispstar * lstar[I]
        pgala[I] = (1. - thispstar) * lgala[I]
        tplot = np.arange(0.8, 10.0, 0.01)
        pplot = probability(tplot, thispstar, scale)
        if q < 3:
            baseline = splits[q]
        else:
            baseline = 25.
        plt.plot(tplot, baseline - 0.7 * pplot / np.max(pplot), 'r-', lw=2., alpha=0.75)
    plt.savefig('%s-data-models.%s' % (prefix, suffix))
    print pstarbest, scalebest

    plt.clf()
    plt.plot(t, m, 'ko', mew=0, alpha=0.5)
    plt.xlim(sizelim)
    plt.xlabel(r'size $\theta$')
    plt.ylim(maglim)
    plt.ylabel('magnitude $m$')
    excut = 1.1
    plt.axvline(excut, color='r', lw=2, alpha=0.75)
    plt.savefig('%s-data-cut.%s' % (prefix, suffix))

    cuts, completeness, purity = straight_cut(t, c)
    plt.clf()
    plt.plot(cuts, completeness, 'k-')
    plt.plot(cuts, purity, 'k--')
    plt.xlim(np.min(cuts), np.max(cuts))
    plt.xlabel(r'size cut $\theta_c$')
    plt.ylim(0., 1.)
    plt.title('hard cut: completeness and purity')
    plt.savefig('%s-hard.%s' % (prefix, suffix))

    plt.clf()
    plt.plot(t, m, 'ko', mew=0, alpha=0.5)
    plt.xlim(sizelim)
    plt.xlabel(r'size $\theta$')
    plt.ylim(maglim)
    plt.ylabel('magnitude $m$')
    plt.plot([excut, excut, -1.], [0., 24., 24.], 'r-', lw=2, alpha=0.75)
    plt.savefig('%s-data-cut-24.%s' % (prefix, suffix))

    cuts, completeness, purity = straight_cut(t[m < 24.], c[m < 24.], ntotal=np.sum(c==0))
    plt.clf()
    plt.plot(cuts, completeness, 'k-')
    plt.plot(cuts, purity, 'k--')
    plt.xlim(np.min(cuts), np.max(cuts))
    plt.xlabel(r'size cut $\theta_c$')
    plt.ylim(0., 1.)
    plt.title('hard cut: completeness and purity for $m < 24$')
    plt.savefig('%s-hard-24.%s' % (prefix, suffix))

    cuts, completeness, purity = prob_cut(np.log(lstar / lgala), c)
    plt.clf()
    plt.plot(cuts, completeness, 'k-')
    plt.plot(cuts, purity, 'k--')
    plt.xlim(np.min(cuts), np.max(cuts))
    plt.xlabel('ln likelihood ratio cut')
    plt.ylim(0., 1.)
    plt.title('likelihood ratio cut: completeness and purity')
    plt.savefig('%s-like.%s' % (prefix, suffix))

    cuts, completeness, purity = prob_cut(np.log(pstar / pgala), c)
    plt.clf()
    plt.plot(cuts, completeness, 'k-')
    plt.plot(cuts, purity, 'k--')
    plt.xlim(np.min(cuts), np.max(cuts))
    plt.xlabel('ln probability ratio cut')
    plt.ylim(0., 1.)
    plt.title('probability ratio cut: completeness and purity')
    plt.savefig('%s-prob.%s' % (prefix, suffix))

if __name__ == '__main__':
    main()

reason: import numpy,import scipy


text:
# -*- coding: utf-8 -*-
"""
Created on Fri Mar 12 19:13:31 2021

@author: aymer
"""

from matplotlib import pyplot as plt
import numpy as np
# Create the data set
x = np.arange(0, 10, 0.05)
y = np.sin(x)
# Define the confidence interval
ci = 0.1 * np.std(y) / np.mean(y)
# Plot the sinus function
plt.plot(x, y)
# Plot the confidence interval
plt.fill_between(x, (y-ci), (y+ci), color='blue', alpha=0.1)
plt.show()

max_stars_repo_name: aymericvie/networked-genetic-algorithms

max_issues_repo_name: aymericvie/networked-genetic-algorithms

max_forks_repo_name: aymericvie/networked-genetic-algorithms


text:
import torch
from torch import nn
from torch import distributions
import numpy as np

import pfrl

from pfrl import explorers
from pfrl.replay_buffer import high_level_batch_experiences_with_goal
from pfrl.agents import HIROHighLevelGoalConditionedTD3, GoalConditionedTD3
from pfrl.nn import ConstantsMult
from pfrl.nn.lmbda import Lambda

class HRLControllerBase():
    def __init__(
            self,
            state_dim,
            goal_dim,
            action_dim,
            scale,
            replay_buffer,
            actor_lr,
            critic_lr,
            expl_noise,
            policy_noise,
            noise_clip,
            gamma,
            policy_freq,
            tau,
            is_low_level,
            buffer_freq,
            minibatch_size,
            gpu,
            add_entropy,
            burnin_action_func=None,
            replay_start_size=2500,
            temperature=1.0,
            optimize_temp=False):
        self.scale = scale

        if gpu is not None and gpu >= 0:
            assert torch.cuda.is_available()
            self.device = torch.device("cuda:{}".format(gpu))
        else:
            self.device = torch.device("cpu")

        self.scale_tensor = torch.tensor(self.scale).float().to(self.device)
        # parameters
        self.expl_noise = expl_noise
        self.policy_noise = policy_noise
        self.noise_clip = noise_clip
        self.gamma = gamma
        self.policy_freq = policy_freq
        self.tau = tau
        self.is_low_level = is_low_level
        self.minibatch_size = minibatch_size
        self.add_entropy = add_entropy

        # create agent
        if self.add_entropy:
            def squashed_diagonal_gaussian_head(x):
                """
                taken from the SAC code.
                """
                assert x.shape[-1] == action_dim * 2

                mean, log_scale = torch.chunk(x, 2, dim=-1)
                log_scale = torch.clamp(log_scale, -20.0, 2.0)
                var = torch.exp(log_scale * 2)
                base_distribution = distributions.Independent(
                    distributions.Normal(loc=mean, scale=torch.sqrt(var)), 1
                )
                # cache_size=1 is required for numerical stability
                return distributions.transformed_distribution.TransformedDistribution(
                    base_distribution, [distributions.transforms.TanhTransform(cache_size=1)]
                )

            # SAC policy definition:
            policy = nn.Sequential(
                nn.Linear(state_dim + goal_dim, 256),
                nn.ReLU(),
                nn.Linear(256, 256),
                nn.ReLU(),
                nn.Linear(256, action_dim * 2),
                Lambda(squashed_diagonal_gaussian_head),
                )

            torch.nn.init.xavier_uniform_(policy[0].weight)
            torch.nn.init.xavier_uniform_(policy[2].weight)
            torch.nn.init.xavier_uniform_(policy[4].weight)
            explorer = explorers.AdditiveGaussian(
                scale=0.0,
            )

        else:
            policy = nn.Sequential(
                nn.Linear(state_dim + goal_dim, 300),
                nn.ReLU(),
                nn.Linear(300, 300),
                nn.ReLU(),
                nn.Linear(300, action_dim),
                nn.Tanh(),
                pfrl.policies.DeterministicHead(),
                )
            # TODO - have proper low and high values from action space.
            # from the hiro paper, the scale is 1.0
            explorer = explorers.AdditiveGaussian(
                scale=self.expl_noise,
                low=-self.scale,
                high=self.scale
            )

        policy_optimizer = torch.optim.Adam(policy.parameters(), lr=actor_lr)

        def make_q_func_with_optimizer():
            q_func = nn.Sequential(
                pfrl.nn.ConcatObsAndAction(),
                nn.Linear(state_dim + goal_dim + action_dim, 300),
                nn.ReLU(),
                nn.Linear(300, 300),
                nn.ReLU(),
                nn.Linear(300, 1),
            )
            q_func_optimizer = torch.optim.Adam(q_func.parameters(), lr=critic_lr)
            return q_func, q_func_optimizer

        q_func1, q_func1_optimizer = make_q_func_with_optimizer()
        q_func2, q_func2_optimizer = make_q_func_with_optimizer()

        def default_target_policy_smoothing_func(batch_action):
            """Add noises to actions for target policy smoothing."""
            noise = torch.clamp(self.policy_noise * torch.randn_like(batch_action), -self.noise_clip, self.noise_clip)
            smoothed_action = batch_action + noise
            smoothed_action = torch.min(smoothed_action, torch.tensor(self.scale).to(self.device).float())
            smoothed_action = torch.max(smoothed_action, torch.tensor(-self.scale).to(self.device).float())
            return smoothed_action

        input_scale = self.scale_tensor

        if self.is_low_level:
            # standard goal conditioned td3
            self.agent = GoalConditionedTD3(
                policy,
                q_func1,
                q_func2,
                policy_optimizer,
                q_func1_optimizer,
                q_func2_optimizer,
                replay_buffer,
                gamma=gamma,
                soft_update_tau=tau,
                explorer=explorer,
                update_interval=1,
                policy_update_delay=policy_freq,
                replay_start_size=replay_start_size,
                buffer_freq=buffer_freq,
                minibatch_size=minibatch_size,
                gpu=gpu,
                add_entropy=self.add_entropy,
                scale=input_scale,
                burnin_action_func=burnin_action_func,
                target_policy_smoothing_func=default_target_policy_smoothing_func,
                entropy_temperature=temperature,
                optimize_temp=optimize_temp
                )
        else:
            self.agent = HIROHighLevelGoalConditionedTD3(
                policy,
                q_func1,
                q_func2,
                policy_optimizer,
                q_func1_optimizer,
                q_func2_optimizer,
                replay_buffer,
                gamma=gamma,
                soft_update_tau=tau,
                explorer=explorer,
                update_interval=1,
                policy_update_delay=policy_freq,
                replay_start_size=replay_start_size/buffer_freq - 5,
                buffer_freq=buffer_freq,
                minibatch_size=minibatch_size,
                gpu=gpu,
                add_entropy=self.add_entropy,
                scale=input_scale,
                burnin_action_func=burnin_action_func,
                target_policy_smoothing_func=default_target_policy_smoothing_func,
                entropy_temperature=temperature,
                optimize_temp=optimize_temp
                )

        self.device = self.agent.device

    def save(self, directory):
        """
        save the internal state of the TD3 agent.
        """
        self.agent.save(directory)

    def load(self, directory):
        """
        load the internal state of the TD3 agent.
        """
        self.agent.load(directory)

    def policy(self, state, goal):
        """
        run the policy (actor).
        """
        action = self.agent.act_with_goal(torch.FloatTensor(state), torch.FloatTensor(goal))
        return np.clip(action, a_min=-self.scale, a_max=self.scale)

    def _observe(self, states, goals, rewards, done, state_arr=None, action_arr=None):
        """
        observe, and train (if we can sample from the replay buffer)
        """
        self.agent.observe_with_goal(torch.FloatTensor(states), torch.FloatTensor(goals), rewards, done, None)

    def observe(self, states, goals, rewards, done, iterations=1):
        """
        get data from the replay buffer, and train.
        """
        return self._observe(states, goals, rewards, goals, done)

# lower controller
class LowerController(HRLControllerBase):
    def __init__(
            self,
            state_dim,
            goal_dim,
            action_dim,
            scale,
            replay_buffer,
            add_entropy,
            actor_lr=0.0001,
            critic_lr=0.001,
            expl_noise=1.0,
            policy_noise=0.2,
            noise_clip=0.5,
            gamma=0.99,
            policy_freq=2,
            tau=0.005,
            is_low_level=True,
            buffer_freq=10,
            minibatch_size=100,
            gpu=None,
            burnin_action_func=None,
            temperature=1.0,
            optimize_low_temp=False):
        super(LowerController, self).__init__(
                                            state_dim=state_dim,
                                            goal_dim=goal_dim,
                                            action_dim=action_dim,
                                            scale=scale,
                                            replay_buffer=replay_buffer,
                                            actor_lr=actor_lr,
                                            critic_lr=critic_lr,
                                            expl_noise=expl_noise,
                                            policy_noise=policy_noise,
                                            noise_clip=noise_clip,
                                            gamma=gamma,
                                            policy_freq=policy_freq,
                                            tau=tau,
                                            is_low_level=is_low_level,
                                            buffer_freq=buffer_freq,
                                            minibatch_size=minibatch_size,
                                            gpu=gpu,
                                            add_entropy=add_entropy,
                                            burnin_action_func=burnin_action_func,
                                            temperature=temperature,
                                            optimize_temp=optimize_low_temp)

    def observe(self, n_s, g, r, done):

        return self._observe(n_s, g, r, done)

# higher controller

class HigherController(HRLControllerBase):
    def __init__(
            self,
            state_dim,
            goal_dim,
            action_dim,
            scale,
            replay_buffer,
            add_entropy,
            actor_lr=0.0001,
            critic_lr=0.001,
            expl_noise=1.0,
            policy_noise=0.2,
            noise_clip=0.5,
            gamma=0.99,
            policy_freq=2,
            tau=0.005,
            is_low_level=False,
            buffer_freq=10,
            minibatch_size=100,
            gpu=None,
            burnin_action_func=None,
            temperature=1.0,
            optimize_high_temp=False):
        super(HigherController, self).__init__(
                                                state_dim=state_dim,
                                                goal_dim=goal_dim,
                                                action_dim=action_dim,
                                                scale=scale,
                                                replay_buffer=replay_buffer,
                                                actor_lr=actor_lr,
                                                critic_lr=critic_lr,
                                                expl_noise=expl_noise,
                                                policy_noise=policy_noise,
                                                noise_clip=noise_clip,
                                                gamma=gamma,
                                                policy_freq=policy_freq,
                                                tau=tau,
                                                is_low_level=is_low_level,
                                                buffer_freq=buffer_freq,
                                                minibatch_size=minibatch_size,
                                                gpu=gpu,
                                                add_entropy=add_entropy,
                                                burnin_action_func=burnin_action_func,
                                                temperature=temperature,
                                                optimize_temp=optimize_high_temp)
        self.action_dim = action_dim

    def _off_policy_corrections(self,
        low_con,
        batch_size,
        states,
        actions,
        next_states,
        low_states,
        low_actions,
        candidate_goals=8):
        # Scale
        goal_dim = self.action_dim
        spec_range = low_con.scale
        # Sample from normal distribution
        loc = (next_states - states)[:, np.newaxis, :goal_dim]
        scale = 0.5 * self.scale[None, None, :]
        original_goal = np.array(actions[:, np.newaxis, :])
        random_goals = np.random.normal(loc=loc, scale=scale, size=(batch_size, candidate_goals, original_goal.shape[-1]))

        candidates = np.concatenate([original_goal, loc, random_goals], axis=1)
        candidates = candidates.clip(-self.scale, self.scale)

        # For ease
        low_actions = np.array(low_actions)
        seq_len = len(low_states[0])
        new_batch_sz = seq_len * batch_size
        low_action_dim = low_actions[0][0].shape
        low_obs_dim = low_states[0][0].shape
        ncands = candidates.shape[1]

        true_low_actions = low_actions.reshape((new_batch_sz,) + low_action_dim)
        observations = low_states.reshape((new_batch_sz,) + low_obs_dim)
        goal_shape = (new_batch_sz, self.action_dim)

        pred_actions = np.zeros((ncands, new_batch_sz) + low_action_dim)

        low_con.agent.training = False
        for c in range(ncands):
            subgoal = candidates[:, c]
            candidate = (subgoal + low_states[:, 0, :self.action_dim])[:, None] - low_states[:, :, :self.action_dim]
            candidate = candidate.reshape(*goal_shape)
            pred_actions[c] = low_con.policy(torch.tensor(observations).float(), torch.tensor(candidate).float())
        low_con.agent.training = True

        difference = (pred_actions - true_low_actions)
        # difference = np.where(difference != -np.inf, difference, 0)
        difference = difference.reshape((ncands, batch_size, seq_len) + low_action_dim).transpose(1, 0, 2, 3)

        normalized_error = - np.square(difference) / np.square(spec_range)
        fitness = np.sum(normalized_error, axis=(2, 3))
        best_actions = np.argmax(fitness, axis=-1)

        return candidates[np.arange(batch_size), best_actions]

    def update(self, low_con):
        batch = self.agent.sample_if_possible()
        if batch:
            experience = high_level_batch_experiences_with_goal(batch, self.device,
                lambda x: x, self.gamma)
            states = experience['state']
            actions = experience['action']
            next_states = experience['next_state']
            action_arr = experience['action_arr']
            state_arr = experience['state_arr']

            actions = self._off_policy_corrections(
                low_con,
                self.minibatch_size,
                states.cpu().data.numpy(),
                actions.cpu().data.numpy(),
                next_states.cpu().data.numpy(),
                state_arr.cpu().data.numpy(),
                action_arr.cpu().data.numpy())

            tensor_actions = torch.FloatTensor(actions).to(self.agent.device)
            # relabel actions
            experience['action'] = tensor_actions

            self.agent.high_level_update_batch(experience)

    def observe(self, state_arr, action_arr, r, g, n_s, done):
        """
        train the high level controller with
        the novel off policy correction.
        """
        # step 1 - record experience in replay buffer
        self.agent.observe_with_goal_state_action_arr(torch.FloatTensor(state_arr),
                                                      torch.FloatTensor(action_arr),
                                                      torch.FloatTensor(n_s),
                                                      torch.FloatTensor(g), r, done, None)

max_stars_repo_name: HIRO-group/pfrl

max_issues_repo_name: HIRO-group/pfrl

max_forks_repo_name: HIRO-group/pfrl


text:
"""
Fast method to get the proper indexes for a circular buffer.

See history of version control. Python C tuple '_circular_indexes.c' building a tuple was much slower.
"""
import numpy as np

__all__ = ['get_indexes']

def get_indexes(start, length, maxsize):
    """Return the indexes from the given start position to the given length."""
    stop = start + length

    if stop > maxsize:
        # Check roll-over/roll-under
        try:
            return np.concatenate((np.arange(start, maxsize),
                                   np.arange(0, stop % maxsize)))
        except ZeroDivisionError:
            return []
    elif stop < 0:
        # Negative length roll-under
        return np.concatenate((np.arange(start, -1, -1),
                               np.arange(maxsize, maxsize-stop, -1)))

    # Return a simple slice
    try:
        return slice(start, stop, length//abs(length))
    except ZeroDivisionError:
        return slice(start, stop)
# end get_indexes


text:
"""
This file is part of the OpenProtein project.

For license information, please see the LICENSE file in the root directory.
"""

import json
import time
import numpy as np
import requests
import torch.optim as optim
from util import set_experiment_id, write_out, write_model_to_disk, write_result_summary

def train_model(data_set_identifier, model, train_loader, validation_loader,
                learning_rate, minibatch_size=64, eval_interval=50, hide_ui=False,
                use_gpu=False, minimum_updates=1000):
    set_experiment_id(data_set_identifier, learning_rate, minibatch_size)

    validation_dataset_size = validation_loader.dataset.__len__()

    if use_gpu:
        model = model.cuda()

    optimizer = optim.Adam(model.parameters(), lr=learning_rate)

    sample_num = list()
    train_loss_values = list()
    validation_loss_values = list()

    best_model_loss = 1e20
    best_model_minibatch_time = None
    best_model_path = None
    best_json_data = None
    stopping_condition_met = False
    minibatches_proccesed = 0

    while not stopping_condition_met:
        optimizer.zero_grad()
        model.zero_grad()
        loss_tracker = np.zeros(0)
        for _minibatch_id, training_minibatch in enumerate(train_loader, 0):
            minibatches_proccesed += 1
            start_compute_loss = time.time()
            loss = model.compute_loss(training_minibatch)
            write_out("Train loss:", float(loss))
            start_compute_grad = time.time()
            loss.backward()
            loss_tracker = np.append(loss_tracker, float(loss))
            end = time.time()
            write_out("Loss time:", start_compute_grad - start_compute_loss, "Grad time:",
                      end - start_compute_grad)
            optimizer.step()
            optimizer.zero_grad()
            model.zero_grad()

            # for every eval_interval samples, plot performance on the validation set
            if minibatches_proccesed % eval_interval == 0:

                write_out("Testing model on validation set...")

                train_loss = float(loss_tracker.mean())
                loss_tracker = np.zeros(0)
                validation_loss, json_data, _ = model.evaluate_model(validation_loader)

                if validation_loss < best_model_loss:
                    best_model_loss = validation_loss
                    best_model_minibatch_time = minibatches_proccesed
                    best_model_path = write_model_to_disk(model)
                    best_json_data = json_data

                write_out("Validation loss:", validation_loss, "Train loss:", train_loss)
                write_out("Best model so far (validation loss): ", best_model_loss, "at time",
                          best_model_minibatch_time)
                write_out("Best model stored at " + best_model_path)
                write_out("Minibatches processed:", minibatches_proccesed)
                sample_num.append(minibatches_proccesed)
                train_loss_values.append(train_loss)
                validation_loss_values.append(validation_loss)

                json_data["validation_dataset_size"] = validation_dataset_size
                json_data["sample_num"] = sample_num
                json_data["train_loss_values"] = train_loss_values
                json_data["validation_loss_values"] = validation_loss_values

                write_out(json_data)

                if not hide_ui:
                    res = requests.post('http://localhost:5000/graph', json=json_data)
                    if res.ok:
                        print(res.json())

                if minibatches_proccesed > minimum_updates and minibatches_proccesed \
                        >= best_model_minibatch_time + minimum_updates:
                    stopping_condition_met = True
                    break
    write_result_summary(best_model_loss)
    write_result_summary(json.dumps(best_json_data))
    return best_model_path

max_stars_repo_name: biokvantz/openprotein

max_issues_repo_name: biokvantz/openprotein

max_forks_repo_name: biokvantz/openprotein