text:
#!/usr/bin/env python

'''
Optimizations over an expression graph.

Optimization passes take as input an expression graph, and return a
(hopefully) simpler, equivalent graph.  This module defines the
pass infrastructure, the fusion passes and an optimization pass to
lower code to Parakeet.
'''
import operator
import math
import weakref

from collections import namedtuple

from . import local
from . import tiling
from .base import Expr, Val, AsArray, ListExpr, lazify, expr_like, ExprTrace
from .base import NotShapeable, CollectionExpr
from .filter import FilterExpr
from .local import LocalInput, LocalMapExpr, LocalMapLocationExpr, make_var
from .local import ParakeetExpr
from .map import MapExpr, Map2Expr
from .ndarray import NdArrayExpr
from .outer import OuterProductExpr
from .shuffle import ShuffleExpr
from .slice import SliceExpr
from .write_array import WriteArrayExpr
from .reduce import ReduceExpr, LocalReduceExpr
from ..dot import DotExpr

from ... import util
from ...util import Assert
from ...config import FLAGS, BoolFlag
from ...array.distarray import DistArray, LocalWrapper

try:
  import numexpr
except:
  numexpr = None

try:
  import parakeet
except:
  parakeet = None

_parakeet_blacklist = set()
_tiled_exprlist = {}
_not_idempotent_list = set()

def disable_parakeet(fn):
  "Disables parakeet optimization for this function."
  _parakeet_blacklist.add(fn)
  return fn

def not_idempotent(fn):
  "Force the result of not_idempotent ``fn`` to be evaluated at once."
  def wrapped(*args, **kw):
    result = fn(*args, **kw)
    if isinstance(result, Expr):
      result.needs_cache = True
      _not_idempotent_list.add(id(result))
    return result
  return wrapped

visited_expr = {'map_fusion': weakref.WeakValueDictionary(),
                'reduce_fusion': weakref.WeakValueDictionary(),
                'collapse_cached': weakref.WeakValueDictionary(),
                'parakeet_gen': weakref.WeakValueDictionary(),
                'rotate_slice': weakref.WeakValueDictionary(),
                'auto_tiling': weakref.WeakValueDictionary()
                }

class OptimizePass(object):
  def __init__(self):
    self.visited = visited_expr[self.name]
    #self.visited = {}

  def visit(self, op):
    if not isinstance(op, Expr):
      return op

    if op.expr_id in self.visited:
      return self.visited[op.expr_id]

    #assert not op in self.visited, 'Infinite recursion during optimization %s' % op
    #self.visited.add(op)

    util.log_debug('VISIT %s: %s', op.typename(), hash(op))
    opt_op = None
    if hasattr(self, 'visit_default'):
      opt_op = self.visit_default(op)
    elif hasattr(self, 'visit_%s' % op.typename()):
      opt_op = getattr(self, 'visit_%s' % op.typename())(op)
    else:
      opt_op = op.visit(self)

    self.visited[opt_op.expr_id] = opt_op
    return opt_op

def fusable(v):
  return isinstance(v, (MapExpr,
                        ReduceExpr,
                        ShuffleExpr,
                        NdArrayExpr,
                        SliceExpr,
                        FilterExpr,
                        Val,
                        AsArray,
                        WriteArrayExpr))

def merge_var(children, child_to_var, k, v):
  """Add a new expression with key ``k`` to the ``children`` dictionary.

  If ``k`` is already in the dictionary, than ``v`` must be equal to
  the current value stored there.
  """
  try:
    i = child_to_var.index(k)
    Assert.eq(v, children[i])
  except ValueError:
    children.append(v)
    child_to_var.append(k)

class MapMapFusion(OptimizePass):
  '''Fold sequences of Map operations together.

  map(f, map(g, map(h, x))) -> map(f . g . h, x)
  '''
  name = 'map_fusion'

  def visit_MapExpr(self, expr):
    Assert.iterable(expr.children)
    #util.log_info('VISIT %d', id(op))
    #util.log_info('VISIT %s', op.children)
    map_children = self.visit(expr.children)
    all_maps = True
    Assert.isinstance(map_children, ListExpr)
    for k, v in zip(expr.child_to_var, map_children):
      if not fusable(v):
        util.log_debug('Skipping fusion: (%s -> %s)' % (k, type(v)))
        all_maps = False
        break

    if not all_maps or isinstance(expr.op, local.ParakeetExpr) or \
       id(expr) in _not_idempotent_list:
      return expr.visit(self)

    #util.log_info('Original: %s', expr.op)
    children = []
    child_to_var = []
    combined_op = expr.op.__class__(fn=expr.op.fn, kw=expr.op.kw,
                                    pretty_fn=expr.op.pretty_fn)

    trace = ExprTrace()

    for child_expr in map_children:
      trace.fuse(child_expr.stack_trace)

      if isinstance(child_expr, MapExpr):
        for k, v in zip(child_expr.child_to_var, child_expr.children):
          merge_var(children, child_to_var, k, v)

        combined_op.add_dep(child_expr.op)
        util.log_debug('Fusion: %s <- %s', expr.expr_id, child_expr.expr_id)
      else:
        children.append(child_expr)
        key = make_var()
        combined_op.add_dep(LocalInput(idx=key))
        child_to_var.append(key)

    if isinstance(combined_op, LocalMapLocationExpr):
      combined_op.add_dep(LocalInput(idx='extent'))

    return expr_like(expr,
                     children=ListExpr(vals=children),
                     child_to_var=child_to_var,
                     op=combined_op,
                     trace=trace)

class ReduceMapFusion(OptimizePass):
  '''Fuse reduce(f, map(g, X)) -> reduce(f . g, X)'''
  name = 'reduce_fusion'

  def visit_ReduceExpr(self, expr):
    Assert.isinstance(expr.children, ListExpr)
    old_children = self.visit(expr.children)

    for v in old_children:
      if not isinstance(v, (MapExpr, ParakeetExpr)) or id(v) in _not_idempotent_list:
        return expr.visit(self)

    combined_op = LocalReduceExpr(fn=expr.op.fn,
                                  kw=expr.op.kw,
                                  deps=[expr.op.deps[0]])

    new_children = []
    new_child_to_var = []
    trace = ExprTrace()
    for i in range(len(old_children)):
      name = expr.child_to_var[i]
      child_expr = old_children[i]
      for j in range(len(child_expr.children)):
        k = child_expr.child_to_var[j]
        v = child_expr.children[j]
        merge_var(new_children, new_child_to_var, k, v)
      combined_op.add_dep(child_expr.op)
      trace.fuse(child_expr.stack_trace)

    return expr_like(expr,
                     children=ListExpr(vals=new_children),
                     child_to_var=new_child_to_var,
                     axis=expr.axis,
                     dtype_fn=expr.dtype_fn,
                     accumulate_fn=expr.accumulate_fn,
                     op=combined_op,
                     tile_hint=expr.tile_hint,
                     trace=trace)

class CollapsedCachedExpressions(OptimizePass):
  '''Replace expressions which have already been evaluated
  with a simple value expression.

  This results in simpler local expressions when evaluating
  iterative programs.
  '''

  name = 'collapse_cached'

  def visit_default(self, expr):
    #util.log_info('Visit: %s, %s', expr.expr_id, expr.cache)
    cache = expr.cache()
    if cache is not None:
      util.log_info('Collapsing %s %s', expr.expr_id, expr.typename())
      return lazify(cache)
    else:
      return expr.visit(self)

def _find_modules(op):
  '''Find any modules referenced by the given `LocalOp` or its dependencies'''
  modules = set()
  if isinstance(op, local.FnCallExpr):
    if hasattr(op.fn, '__module__'):
      modules.add(op.fn.__module__)
    elif hasattr(op.fn, '__class__'):
      modules.add(op.fn.__class__.__module__)

  for d in op.deps:
    modules.union(_find_modules(d))

  return modules

def _parakeet_codegen(op):
  '''Given a local operation, generate an equivalent parakeet function definition.'''

  def _codegen(op):
    if isinstance(op, local.FnCallExpr):
      if util.is_lambda(op.fn) and not op.pretty_fn:
        raise local.CodegenException('Cannot codegen through a lambda expression: %s' % op.fn)

      if op.fn in _parakeet_blacklist:
        raise local.CodegenException('Blacklisted %s', op.fn)

      name = op.fn_name()

      arg_str = ','.join([_codegen(v) for v in op.deps])
      kw_str = ','.join(['%s=%s' % (k, repr(v)) for k, v in op.kw.iteritems()])
      if arg_str:
        kw_str = ',' + kw_str

      return '%s(%s %s)' % (name, arg_str, kw_str)
    elif isinstance(op, local.LocalInput):
      return op.idx
    else:
      raise local.CodegenException('Cannot codegen for %s' % type(op))

  if isinstance(op, ParakeetExpr):
    return op.source

  op_code = _codegen(op)

  module_prelude = ['import parakeet',
                    'import spartan.expr',
                    'import numpy',
                    'from spartan.expr import mathlib',
                    'from spartan import util',
                    ]

  for mod in _find_modules(op):
    module_prelude.append('import %s' % mod)

  fn_prelude = '''
@util.synchronized
@parakeet.jit
'''

  fn = '\n'.join(module_prelude)
  fn = fn + fn_prelude
  fn = fn + 'def _jit_fn'
  fn = fn + '(%s):\n  ' % ','.join(op.input_names())
  fn = fn + 'return ' + op_code
  fn = fn + '\n\n'

  # verify we can compile before proceeding
  local.compile_parakeet_source(fn)
  return fn

class ParakeetGeneration(OptimizePass):
  '''
  Replace local map/reduce operations with an equivalent
  parakeet function definition.
  '''

  name = 'parakeet_gen'
  after = [MapMapFusion, ReduceMapFusion]

  def _should_run_parakeet(self, expr):
    has_fncallexpr = False
    for dep in expr.op.deps:
      if isinstance(dep, local.FnCallExpr):
        has_fncallexpr = True
        break

    # If all deps are not FnCallExprs and len(deps) is a small number(<=2 for now)
    # it is not worth to do code generation.
    if not has_fncallexpr and len(expr.op.deps) <= 2:
      return False
    else:
      return True

  def visit_MapExpr(self, expr):
    # if we've already converted this to parakeet, stop now
    if isinstance(expr.op, local.ParakeetExpr):
      return expr.visit(self)

    if not self._should_run_parakeet(expr):
      return expr.visit(self)

    try:
      source = _parakeet_codegen(expr.op)

      # Recursively check if any children nodes can be optimized with Parakeet
      new_children = []
      for child in expr.children:
        new_children.append(self.visit(child))

      parakeet_expr = expr_like(expr,
                                op=local.ParakeetExpr(source=source, deps=expr.op.deps),
                                children=ListExpr(vals=new_children),
                                child_to_var=expr.child_to_var)

      if id(expr) in _not_idempotent_list: _not_idempotent_list.add(id(parakeet_expr))

      return parakeet_expr
    except local.CodegenException:
      util.log_info('Failed to convert to parakeet.')
      return expr.visit(self)

#   Parakeet doesn't support taking the current Extent object as
#   an argument, which prevents enabling it for reductions.
#
#   def visit_ReduceExpr(self, expr):
#     # if we've already converted this to parakeet, stop now
#     if isinstance(expr.op, local.ParakeetExpr):
#       return expr.visit(self)
#
#     try:
#       source = codegen(expr.op)
#       return expr_like(expr,
#                        children=expr.children,
#                        axis=expr.axis,
#                        dtype_fn=expr.dtype_fn,
#                        accumulate_fn=expr.accumulate_fn,
#                        op=local.ParakeetExpr(source=source,  deps=expr.op.deps))
#     except local.CodegenException:
#       util.log_info('Failed to convert to parakeet.')
#       return expr.visit(self)

class RotateSlice(OptimizePass):
  '''
  This pass rotates slice operations to the bottom of the expression graph.

  By moving slices down, we open up opportunities for performing more
  fusion operations.

  The actual operation performed is:

  (a + b)[slice] -> broadcast(a, shape)[slice] + broadcast(b, shape)[slice]
  '''
  name = 'rotate_slice'
  rotated = {}

  def visit_SliceExpr(self, slice_expr):
    'Rotate this slice with a child map expression.'
    map_expr = slice_expr.src

    if not isinstance(map_expr, MapExpr):
      return slice_expr.visit(self)

    try:
      map_shape = map_expr.compute_shape()
    except NotShapeable:
      return slice_expr.visit(self)

    Assert.iterable(map_expr.children)
    map_children = self.visit(map_expr.children)

    children = []
    for child_expr in map_children:
      child = SliceExpr(src=child_expr,
                        idx=slice_expr.idx,
                        broadcast_to=map_shape)
      if isinstance(child_expr, MapExpr):
        # If the child is a Map, RotateSlice should traverse it again. Some
        # nodes may be traversed several times. For current implementation,
        # traversing a node several times dones't cause any problem for
        # RotationSlice.
        if self.visited.get(child_expr, None) is not None:
          del self.visited[child_expr]
        child = self.visit(child)

      children.append(child)

    if self.rotated.get(map_expr, None):
      # If this Map has been rotated, we should create a new MapExpr.
      # An example is :
      #    c = a + b
      #    d = c[1:50]
      #    e = c[2:51]
      # In this example, RotateSlice first rotate c[2:51] to a and b.
      # When RotateSlice rotates c[1:50] to a and b, it should create
      # a new Map expr(+) for a[1:50] and b[1:50]
      return MapExpr(children=ListExpr(vals=children),
                     op=map_expr.op,
                     child_to_var=map_expr.child_to_var)
    else:
      self.rotated[map_expr] = True
      return expr_like(map_expr,
                       op=map_expr.op,
                       children=ListExpr(vals=children),
                       child_to_var=map_expr.child_to_var,
                       trace=map_expr.stack_trace)

class AutomaticTiling(OptimizePass):
  '''
  Automatically partition all the arrays.
  We build a min cost max flow DAG for the expr graph, where the cost is the communication cost
  through network. We estimate the cost according to the behavior of each expr. For _builtin_ expr
  and simple map and reduce expr, we can easily estimate the cost. However, for the user defined shuffle
  expr, we can only guess the cost or let users define the cost for us.

  All Exprs:
    [Val, AsArray, DistArray]: Already partitioned array
    [NdArrayExpr]: new array needs to be partitioned
    CollectionExpr([DictExpr, ListExpr, TupleExpr]): self.vals

    [MapExpr]: self.children
    [ReduceExpr]: self.children
    [ShuffleExpr]: self.array, self.fn_kw
    [DotExpr]: self.matrix_a, self.matrix_b

    [SliceExpr, FilterExpr, CheckpointExpr, TileOpExpr]: self.src or self.array
    [WriteArrayExpr]: self.array, self.data

    [TransposeExpr, ReshapeExpr]: self.array
  '''

  name = 'auto_tiling'
  node_type = namedtuple('node_type', ['expr', 'tiling', 'children', 'parents'])
  num_node_per_group = 4
  inited = False
  cost_model = {'map': {(0, 0): 0, (0, 1): 1, (0, 2): 1,
                        (0, 3): 1, (0, 4): 1, (0, 5): 2,
                        (1, 0): 1, (1, 1): 0, (1, 2): 1,
                        (1, 3): 1, (1, 4): 2, (1, 5): 1,
                        (2, 0): 1, (2, 1): 1, (2, 2): 0,
                        (2, 3): 2, (2, 4): 1, (2, 5): 1,
                        (3, 0): 0, (3, 1): 0, (3, 2): 1,
                        (3, 3): 0, (3, 4): 1, (3, 5): 1,
                        (4, 0): 0, (4, 1): 1, (4, 2): 0,
                        (4, 3): 1, (4, 4): 0, (4, 5): 1,
                        (5, 0): 1, (5, 1): 0, (5, 2): 0,
                        (5, 3): 1, (5, 4): 1, (5, 5): 0},
                'map2': {(0, 0): 0, (0, 1): 1, (0, 2): 1, (0, -1): 1,
                         (1, 0): 1, (1, 1): 0, (1, 2): 1, (1, -1): 1,
                         (2, 0): 1, (2, 1): 1, (2, 2): 0, (2, -1): 1,
                         (3, 0): 0, (3, 1): 0, (3, 2): 1, (3, -1): 1,
                         (4, 0): 0, (4, 1): 1, (4, 2): 0, (4, -1): 1,
                         (5, 0): 1, (5, 1): 0, (5, 2): 0, (5, -1): 1}
                }

  def init(self, expr):
    self.cur_node_id = 1
    self.edges = {}
    self.nodes = {0: self.node_type([], -1, [], [])}
    self.expr_to_nodes = {}
    self.split_nodes = {}
    self.groups = []
    self.init_expr = id(expr)
    self.inited = True

    self.tiled_exprlist = _tiled_exprlist
    #self.tiled_exprlist = {}

  def add_edge(self, edge_from, edge_to, edge_cost=0):
    #util.log_warn('add_edge:%d %d cost:%d', edge_from, edge_to, edge_cost)
    if (edge_from, edge_to) not in self.edges:
      self.nodes[edge_from].parents.append(edge_to)
      self.nodes[edge_to].children.append(edge_from)
    self.edges[(edge_from, edge_to)] = edge_cost

  def remove_edge(self, edge_from, edge_to):
    del self.edges[(edge_from, edge_to)]
    self.nodes[edge_from].parents.remove(edge_to)
    self.nodes[edge_to].children.remove(edge_from)

  def add_split_nodes(self, nodes):
    group_id = len(self.groups)
    for node in nodes:
      self.split_nodes[node] = group_id
    self.groups.append(tuple(nodes))
    #util.log_warn('add_split_nodes:%s', nodes)

  def visit_children(self, children, except_child=None):
    child_ids = []
    for child in children:
      if isinstance(child, (Expr, DistArray)) and id(child) != id(except_child):
        child_ids.extend(self.visit_default(child))
    return child_ids

  def visit_NdArrayExpr(self, expr):
    # new array need to be partitioned
    if len(expr.shape) > 1 and expr.shape[1] > 1:
      new_nodes = []

      self.nodes[self.cur_node_id] = self.node_type([expr], 0, [], [])
      new_nodes.append(self.cur_node_id)
      self.add_edge(0, self.cur_node_id, 0)
      self.cur_node_id += 1

      self.nodes[self.cur_node_id] = self.node_type([expr], 1, [], [])
      new_nodes.append(self.cur_node_id)
      self.add_edge(0, self.cur_node_id, 0)
      self.cur_node_id += 1

      self.nodes[self.cur_node_id] = self.node_type([expr], 2, [], [])
      new_nodes.append(self.cur_node_id)
      self.add_edge(0, self.cur_node_id, 0)
      self.cur_node_id += 1

      cost = reduce(operator.mul, expr.shape, 1)
      for i in range(3, self.num_node_per_group):
        self.nodes[self.cur_node_id] = self.node_type([expr], i, [], [])
        new_nodes.append(self.cur_node_id)
        self.add_edge(0, self.cur_node_id, cost)
        self.cur_node_id += 1

      self.add_split_nodes(new_nodes)
      return new_nodes
    else:
      self.nodes[self.cur_node_id] = self.node_type([expr], 0, [], [])
      self.add_edge(0, self.cur_node_id, 0)
      self.cur_node_id += 1
      return [self.cur_node_id - 1]

  def visit_MapExpr(self, expr):
    largest = max(expr.children.vals, key=lambda v: reduce(operator.mul, v.shape, 1))

    child_ids = self.visit_children([largest])
    other_child_ids = self.visit_children(expr.children.vals, largest)
    kw_ids = self.visit_children(expr.op.kw['fn_kw'].itervalues()) if 'fn_kw' in expr.op.kw else []

    # one input map, reuse child expr
    if len(other_child_ids) == 0:
      for child_id in child_ids: self.nodes[child_id].expr.append(expr)
      return child_ids

    if child_ids[0] in self.split_nodes:
      tiling_types = range(self.num_node_per_group)
      self.add_split_nodes(range(self.cur_node_id, self.cur_node_id + self.num_node_per_group))
      expr_node_ids = range(self.cur_node_id, self.cur_node_id + self.num_node_per_group)
    else:
      tiling_types = (self.nodes[child_ids[0]].tiling,)
      expr_node_ids = [self.cur_node_id]

    for (tiling_type, map_child_id) in zip(tiling_types, child_ids):
      self.nodes[self.cur_node_id] = self.node_type([expr], tiling_type, [], [])
      self.add_edge(map_child_id, self.cur_node_id, 0)

      for child_id in other_child_ids:
        child = self.nodes[child_id]
        e_cost = self.cost_model['map'][(child.tiling, tiling_type)] * reduce(operator.mul, child.expr[0].shape, 1)
        self.add_edge(child_id, self.cur_node_id, e_cost)

      for child_id in kw_ids:
        self.add_edge(child_id, self.cur_node_id, reduce(operator.mul, self.nodes[child_id].expr[0].shape, 1))
      self.cur_node_id += 1

    return expr_node_ids

  def visit_ReduceExpr(self, expr):
    child_ids = self.visit_children(expr.children.vals)
    cost = reduce(operator.mul, expr.shape, 1)
    self.nodes[self.cur_node_id] = self.node_type([expr], 0, [], [])
    for child_id in child_ids:
      child = self.nodes[child_id]
      e_cost = 0 if expr.axis is None or child.tiling == 3 or (1-expr.axis) == (child.tiling % 4) else cost
      self.add_edge(child_id, self.cur_node_id, e_cost)
    self.cur_node_id += 1
    return [self.cur_node_id - 1]

  def visit_Map2Expr(self, expr):
    child_id_groups = []
    for array in expr.arrays:
      child_id_groups.append(self.visit_children([array]))

    copy_nodes = []
    for axis, child_ids in zip(expr.axes, child_id_groups):
      if isinstance(axis, tuple): axis = 2
      if axis is None: axis = -1
      self.nodes[self.cur_node_id] = self.node_type([expr], axis, [], [])

      cost = reduce(operator.mul, self.nodes[child_ids[0]].expr[0].shape, 1)
      for child_id in child_ids:
        child = self.nodes[child_id]
        e_cost = self.cost_model['map2'][(child.tiling, axis)] * cost
        self.add_edge(child_id, self.cur_node_id, e_cost)

      if len(copy_nodes) > 0: self.add_edge(self.cur_node_id, copy_nodes[0], 0)
      copy_nodes.append(self.cur_node_id)
      self.cur_node_id += 1

    e_cost = reduce(operator.mul, expr.shape, 1)
    inter_node_id = copy_nodes[0]
    if len(expr.shape) > 1 and expr.shape[1] > 1:
      child_ids = []
      for tiling_type in range(self.num_node_per_group):
        self.nodes[self.cur_node_id] = self.node_type([expr], tiling_type, [], [])
        self.add_edge(inter_node_id, self.cur_node_id, (tiling_type / 3 + 1) * e_cost)
        child_ids.append(self.cur_node_id)
        self.cur_node_id += 1
      self.add_split_nodes(child_ids)
      return child_ids
    else:
      self.nodes[self.cur_node_id] = self.node_type([expr], 0, [], [])
      self.add_edge(inter_node_id, self.cur_node_id, e_cost)
      self.cur_node_id += 1
      return [self.cur_node_id - 1]

  def visit_OuterProductExpr(self, expr):
    child_id_groups = []
    for array in expr.arrays:
      child_id_groups.append(self.visit_children([array]))

    copy_nodes = []
    for axis, child_ids in zip(expr.axes, child_id_groups):
      self.nodes[self.cur_node_id] = self.node_type([expr], axis, [], [])

      cost = reduce(operator.mul, self.nodes[child_ids[0]].expr[0].shape, 1)
      for child_id in child_ids:
        child = self.nodes[child_id]
        e_cost = 0 if axis is not None and (axis == (child.tiling % 4) or child.tiling == 3) else cost
        self.add_edge(child_id, self.cur_node_id, e_cost)

      if len(copy_nodes) > 0: self.add_edge(self.cur_node_id, copy_nodes[0], cost)
      copy_nodes.append(self.cur_node_id)
      self.cur_node_id += 1

    e_cost = reduce(operator.mul, expr.shape, 1)
    inter_node_id = copy_nodes[0]
    if len(expr.shape) > 1 and expr.shape[1] > 1:
      child_ids = []
      for tiling_type in range(self.num_node_per_group):
        self.nodes[self.cur_node_id] = self.node_type([expr], tiling_type, [], [])
        self.add_edge(inter_node_id, self.cur_node_id, (tiling_type / 3 + 1) * e_cost)
        child_ids.append(self.cur_node_id)
        self.cur_node_id += 1
      self.add_split_nodes(child_ids)
      return child_ids
    else:
      self.nodes[self.cur_node_id] = self.node_type([expr], 0, [], [])
      self.add_edge(inter_node_id, self.cur_node_id, e_cost)
      self.cur_node_id += 1
      return [self.cur_node_id - 1]

  def visit_ShuffleExpr(self, expr):
    for child in expr.fn_kw.itervalues():
      if isinstance(child, (Expr, DistArray)) and hash(child) not in expr.cost_hint:
        cost = reduce(operator.mul, child.shape, 1)
        expr.cost_hint[hash(child)] = {'00': cost, '01': cost, '10': cost, '11': cost}
    if expr.target is not None and hash(expr.target) not in expr.cost_hint:
      cost = reduce(operator.mul, expr.target.shape, 1)
      expr.cost_hint[hash(expr.target)] = {'00': cost, '01': cost, '10': cost, '11': cost}

    child_ids = self.visit_children([expr.array])
    other_child_ids = self.visit_children(expr.fn_kw.itervalues())

    # calc the copy cost
    if len(other_child_ids) == 0:
      expr_node_ids = child_ids
    else:
      if child_ids[0] in self.split_nodes:
        tiling_types = (0, 1)
        self.add_split_nodes(self.cur_node_id, self.cur_node_id + 1)
        expr_node_ids = [self.cur_node_id, self.cur_node_id + 1]
      else:
        tiling_types = (self.nodes[child_ids[0]].tiling,)
        expr_node_ids = [self.cur_node_id]

      for i in xrange(len(tiling_types)):
        tiling_type = tiling_types[i]
        self.nodes[self.cur_node_id] = self.node_type([expr], tiling_type, [], [])
        self.add_edge(child_ids[i], self.cur_node_id, 0)

        for child in expr.fn_kw.itervalues():
          if isinstance(child, (Expr, DistArray)):
            for child_id in self.expr_to_nodes[hash(child)]:
              self.add_edge(child_id, self.cur_node_id, expr.cost_hint[hash(child)]['%d%d' % (self.nodes[child_id].tiling, tiling_type)])
        self.cur_node_id += 1

    # calculate update cost
    if expr.target is not None:
      other_child_ids = expr_node_ids
      child_ids = self.visit_children([expr.target])

      if child_ids[0] in self.split_nodes:
        tiling_types = (0, 1)
        self.add_split_nodes(self.cur_node_id, self.cur_node_id + 1)
        expr_node_ids = [self.cur_node_id, self.cur_node_id + 1]
      else:
        tiling_types = (self.nodes[child_ids[0]].tiling,)
        expr_node_ids = [self.cur_node_id]

      for i in xrange(len(tiling_types)):
        tiling_type = tiling_types[i]
        self.nodes[self.cur_node_id] = self.node_type([expr], tiling_type, [], [])
        self.add_edge(child_ids[i], self.cur_node_id, 0)

        for child_id in other_child_ids:
          self.add_edge(child_id, self.cur_node_id, expr.cost_hint[hash(expr.target)]['%d%d' % (self.nodes[child_id].tiling, tiling_type)])
        self.cur_node_id += 1
      return expr_node_ids

    if expr.shape == expr.array.shape or len(other_child_ids) > 0:
      if len(other_child_ids) == 0:
        for child_id in expr_node_ids: self.nodes[child_id].expr.append(expr)
      return expr_node_ids

    # fix result tiling
    if len(expr.shape) <= 1 or 1 in expr.shape[:2]:
      tiling_type = 0 if len(expr.shape) <= 1 else 1 - expr.shape.index(1)
      self.nodes[self.cur_node_id] = self.node_type([expr], tiling_type, [], [])
      for child_id in expr_node_ids: self.add_edge(child_id, self.cur_node_id, 0)
      self.cur_node_id += 1
      return [self.cur_node_id - 1]
    else:
      child_ids = expr_node_ids
      expr_node_ids = []
      for child_id in child_ids:
        self.nodes[self.cur_node_id] = self.node_type([expr], self.nodes[child_id].tiling, [], [])
        expr_node_ids.append(self.cur_node_id)
        self.add_edge(child_id, self.cur_node_id, 0)
        self.cur_node_id += 1
      if len(expr_node_ids) > 1: self.add_split_nodes(*expr_node_ids)
    return expr_node_ids

  def visit_DotExpr(self, expr):
    child_ids = self.visit_children([expr.matrix_a])
    other_child_ids = self.visit_children([expr.matrix_b])

    # calc copy cost
    if len(other_child_ids) != 0:
      if child_ids[0] in self.split_nodes:
        tiling_types = range(self.num_node_per_group)
        self.add_split_nodes(range(self.cur_node_id, self.cur_node_id + self.num_node_per_group))
        expr_node_ids = range(self.cur_node_id, self.cur_node_id + self.num_node_per_group)
      else:
        tiling_types = (self.nodes[child_ids[0]].tiling,)
        expr_node_ids = [self.cur_node_id]

      for i in xrange(len(tiling_types)):
        tiling_type = tiling_types[i]
        self.nodes[self.cur_node_id] = self.node_type([expr], tiling_type, [], [])
        self.add_edge(child_ids[i], self.cur_node_id, 0)

        for child_id in other_child_ids:
          child = self.nodes[child_id]
          e_cost = 0 if child.tiling in (0, 3) and tiling_type in (1, 2) else reduce(operator.mul, child.expr[0].shape, 1)
          self.add_edge(child_id, self.cur_node_id, e_cost)
        self.cur_node_id += 1

      child_ids = expr_node_ids

    # calc update cost
    if len(expr.shape) == 1 or expr.shape[1] == 1:
      tiling_types = (0,)
      expr_node_ids = [self.cur_node_id]
    else:
      tiling_types = range(self.num_node_per_group)
      self.add_split_nodes(range(self.cur_node_id, self.cur_node_id + self.num_node_per_group))
      expr_node_ids = range(self.cur_node_id, self.cur_node_id + self.num_node_per_group)

    cost = reduce(operator.mul, expr.shape, 1)
    for tiling_type in tiling_types:
      self.nodes[self.cur_node_id] = self.node_type([expr], tiling_type, [], [])
      for child_id in child_ids:
        e_cost = 0 if self.nodes[child_id].tiling in (0, 3) and tiling_type == 0 else cost
        self.add_edge(child_id, self.cur_node_id, e_cost)
      self.cur_node_id += 1
    return expr_node_ids

  def visit_WriteArrayExpr(self, expr):
    child_ids = self.visit_children([expr.array])
    if isinstance(expr.data, (Expr, DistArray)):
      data_child_ids = self.visit_children([expr.data])
      if child_ids[0] in self.split_nodes:
        tiling_types = range(self.num_node_per_group)
        self.add_split_nodes(range(self.cur_node_id, self.cur_node_id + self.num_node_per_group))
        expr_node_ids = range(self.cur_node_id, self.cur_node_id + self.num_node_per_group)
      else:
        tiling_types = (self.nodes[child_ids[0]].tiling,)
        expr_node_ids = [self.cur_node_id]

      for (tiling_type, map_child_id) in zip(tiling_types, child_ids):
        self.nodes[self.cur_node_id] = self.node_type([expr], tiling_type, [], [])
        self.add_edge(map_child_id, self.cur_node_id, 0)

        for child_id in data_child_ids:
          child = self.nodes[child_id]
          e_cost = self.cost_model['map'][(child.tiling, tiling_type)] * reduce(operator.mul, child.expr[0].shape, 1)
          self.add_edge(child_id, self.cur_node_id, e_cost)

        self.cur_node_id += 1
      return expr_node_ids
    else:
      for child_id in child_ids: self.nodes[child_id].expr.append(expr)
      return child_ids

  def visit_aligned_nodes(self, expr, reverse_cost=False):
    array = expr.src if hasattr(expr, 'src') else expr.array
    child_ids = self.visit_children([array])
    if child_ids[0] in self.split_nodes:
      tiling_types = range(self.num_node_per_group)
      self.add_split_nodes(range(self.cur_node_id, self.cur_node_id + self.num_node_per_group))
      expr_node_ids = range(self.cur_node_id, self.cur_node_id + self.num_node_per_group)

      for tiling_type in tiling_types:
        self.nodes[self.cur_node_id] = self.node_type([expr], tiling_type, [], [])
        orig_tiling = -1
        if tiling_type in (2, 3):
          orig_tiling = tiling_type
        elif tiling_type < 2:
          orig_tiling = reverse_cost ^ tiling_type
        else:
          orig_tiling = (reverse_cost ^ (tiling_type % 4)) + 4
        self.add_edge(child_ids[orig_tiling], self.cur_node_id, 0)
        self.cur_node_id += 1
    else:
      child_tiling = self.nodes[child_ids[0]].tiling
      if child_tiling in (2, 3):
        tiling_type = child_tiling
      elif child_tiling < 2:
        tiling_type = reverse_cost ^ child_tiling
      else:
        tiling_type = (reverse_cost ^ (child_tiling % 4)) + 4
      self.nodes[self.cur_node_id] = self.node_type([expr], tiling_type, [], [])
      self.add_edge(child_ids[0], self.cur_node_id, 0)
      expr_node_ids = [self.cur_node_id]
      self.cur_node_id += 1

    return expr_node_ids

  def visit_TransposeExpr(self, expr):
    return self.visit_aligned_nodes(expr, reverse_cost=True)

  def visit_ReshapeExpr(self, expr):
    return self.visit_aligned_nodes(expr, reverse_cost=True)

  def visit_SliceExpr(self, expr):
    return self.visit_aligned_nodes(expr)

  def visit_FilterExpr(self, expr):
    child_ids = self.visit_children([expr.src])
    for child_id in child_ids: self.nodes[child_id].expr.append(expr)
    return child_ids

  def visit_CheckpointExpr(self, expr):
    child_ids = self.visit_children([expr.src])
    for child_id in child_ids: self.nodes[child_id].expr.append(expr)
    return child_ids

  def visit_TileOpExpr(self, expr):
    child_ids = self.visit_children([expr.array])
    for child_id in child_ids: self.nodes[child_id].expr.append(expr)
    return child_ids

  def generate_edges(self, s=0):
    edges = []
    self.nodes[s].parents.sort(key=lambda x: reduce(operator.mul, self.nodes[x].expr[0].shape, 1))
    for parent_id in self.nodes[s].parents:
      edges.append((s, parent_id, self.edges[(s, parent_id)]))
      if parent_id not in self.visited_nodes:
        edges.extend(self.generate_edges(parent_id))
        self.visited_nodes.add(parent_id)
    return edges

  def tile_expr(self, expr, tiling):
    if isinstance(expr, (NdArrayExpr, ReduceExpr, Map2Expr, OuterProductExpr)) and len(expr.shape) > 0:
      expr.tile_hint = list(expr.shape)
      if tiling >= 3:  # duplicate tiling
        print 'dup_tiling', tiling
      elif tiling == 2 and len(expr.shape) > 1:  # block tiling
        expr.tile_hint[0] = int(math.ceil(float(expr.tile_hint[0]) / math.sqrt(FLAGS.num_workers)))
        expr.tile_hint[1] = int(math.ceil(float(expr.tile_hint[1]) / math.sqrt(FLAGS.num_workers)))
        print 'block_tiling', expr.tile_hint, expr.expr_id
      elif len(expr.shape) > tiling:
        expr.tile_hint[tiling] = int(math.ceil(float(expr.tile_hint[tiling]) / FLAGS.num_workers))

  def calc_tiling(self, expr):
    # add T node for graph
    self.nodes[self.cur_node_id] = self.node_type([expr], -1, [], [])
    self.add_edge(self.cur_node_id - 1, self.cur_node_id, 0)
    if self.cur_node_id - 1 in self.split_nodes:
      for i in range(2, self.num_node_per_group + 1):
        self.add_edge(self.cur_node_id - i, self.cur_node_id, 0)
    self.cur_node_id += 1

    # compute best tiling for all exprs
    self.visited_nodes = set()
    edges = self.generate_edges()
    util.log_debug('num of groups %d' % len(self.groups))

    nodes = []
    if FLAGS.tiling_alg == 'maxedge':
      nodes = tiling.maxedge_tiling(self.cur_node_id - 1, edges, self.groups)
      #print 'maxedge', nodes
    elif FLAGS.tiling_alg == 'mincost':
      nodes = tiling.mincost_tiling(self.cur_node_id - 1, edges, self.groups)
      #print 'mincost', nodes
    elif FLAGS.tiling_alg == 'best':
      nodes = tiling.best_tiling(self.cur_node_id - 1, edges, self.groups)
      #print 'best', nodes
    elif FLAGS.tiling_alg == 'worse':
      nodes = tiling.worse_tiling(self.cur_node_id - 1, edges, self.groups)
      #print 'worse', nodes
    else:
      assert 'Unknown tiling algorithm'

    # give expr the best tiling hint
    for node_id in nodes:
      node = self.nodes[node_id]
      if node.tiling < 0 or node.tiling is None: continue
      for cur_expr in node.expr:
        self.tiled_exprlist[hash(cur_expr)] = node.tiling
        self.tile_expr(cur_expr, node.tiling)

    self.inited = False
    return expr

  def tile_cached_expr(self, expr):
    if not isinstance(expr, Expr) or isinstance(expr, (Val, AsArray)): return

    self.tile_expr(expr, self.tiled_exprlist[hash(expr)])

    if hasattr(expr, 'array'): self.tile_cached_expr(expr.array)
    if hasattr(expr, 'src'): self.tile_cached_expr(expr.src)
    if hasattr(expr, 'target'): self.tile_cached_expr(expr.target)
    if isinstance(expr, DotExpr):
      self.tile_cached_expr(expr.matrix_a)
      self.tile_cached_expr(expr.matrix_b)

    if hasattr(expr, 'children'):
      for child in expr.children.vals:
        self.tile_cached_expr(child)
    if hasattr(expr, 'fn_kw') and expr.fn_kw is not None:
      for child in expr.fn_kw.itervalues():
        self.tile_cached_expr(child)
    if hasattr(expr, 'op') and 'fn_kw' in expr.op.kw:
      for child in expr.op.kw['fn_kw'].itervalues():
        self.tile_cached_expr(child)

  def visit_default(self, expr):
    if not self.inited: self.init(expr)
    if hash(expr) in self.tiled_exprlist:
      self.tile_cached_expr(expr)
      tiling = self.tiled_exprlist[hash(expr)]
      self.nodes[self.cur_node_id] = self.node_type([expr], tiling, [], [])
      expr_node_ids = [self.cur_node_id]
      self.add_edge(0, self.cur_node_id, 0)
      self.cur_node_id += 1

    elif hash(expr) in self.expr_to_nodes:
      # cached expr
      expr_node_ids = self.expr_to_nodes[hash(expr)]
      for node_id in expr_node_ids:
        node = self.nodes[node_id]
        node.expr.append(expr)

    elif isinstance(expr, DistArray) or (isinstance(expr, (Val, AsArray)) and isinstance(expr.val, DistArray)):
      # already partitioned array
      array = expr if isinstance(expr, DistArray) else expr.val

      if isinstance(array, LocalWrapper):
        expr_node_ids = []
      else:
        tile_shape = array.tile_shape()
        tiling = 2
        for i in range(len(tile_shape)):
          if tile_shape[i] == array.shape[i]:
            tiling = 1 - i
            break

        self.nodes[self.cur_node_id] = self.node_type([expr], tiling, [], [])
        expr_node_ids = [self.cur_node_id]
        self.add_edge(0, self.cur_node_id, 0)
        self.cur_node_id += 1

    elif isinstance(expr, CollectionExpr):
      # DictExpr, ListExpr, TupleExpr
      children = expr.itervalues() if expr.typename() == 'DictExpr' else expr.vals
      child_ids = self.visit_children(children)
      self.nodes[self.cur_node_id] = self.node_type([expr], -1, [], [])
      expr_node_ids = [self.cur_node_id]
      for child_id in child_ids:
        self.add_edge(child_id, self.cur_node_id, 0)
      self.cur_node_id += 1

    elif hasattr(self, 'visit_%s' % expr.typename()):
      expr_node_ids = getattr(self, 'visit_%s' % expr.typename())(expr)

    else:
      util.log_debug("Skip expr:%s", expr.typename())
      expr_node_ids = []

    # add T node for graph and compute the min cost flow
    if id(expr) == self.init_expr: return self.calc_tiling(expr)

    self.expr_to_nodes[hash(expr)] = expr_node_ids
    return expr_node_ids

def apply_pass(klass, dag):
  if not getattr(FLAGS, 'opt_' + klass.name):
    util.log_debug('Pass %s disabled', klass.name)
    return dag

  util.log_debug('Starting pass %s', klass.name)
  p = klass()
  result = p.visit(dag)
  util.log_debug('Finished pass %s', klass.name)
  return result

passes = []

def optimize(dag):
  if not FLAGS.optimization:
    util.log_info('Optimizations disabled')
    return dag

  util.log_info('Optimization: applying %d passes', len(passes))
  for p in passes:
    dag = apply_pass(p, dag)

  return dag

def add_optimization(klass, default):
  passes.append(klass)

  flagname = 'opt_' + klass.name
  #setattr(Flags, flagname, add_bool_flag(flagname, default=default))
  FLAGS.add(BoolFlag(flagname, default=default, help='Enable %s optimization' % klass.__name__))

  #util.log_info('Passes: %s', passes)

add_optimization(CollapsedCachedExpressions, True)
add_optimization(AutomaticTiling, True)
add_optimization(RotateSlice, False)
add_optimization(MapMapFusion, True)
if parakeet is not None:
  add_optimization(ParakeetGeneration, False)
add_optimization(ReduceMapFusion, True)

FLAGS.add(BoolFlag('optimization', default=True))

max_stars_repo_name: spartan-array/spartan

max_issues_repo_name: spartan-array/spartan

max_forks_repo_name: spartan-array/spartan

reason: import numpy,import numexpr


text:
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import numpy as np
import gc
from datetime import datetime

import matplotlib.pyplot as plt
from pylab import meshgrid,cm,imshow,contour,clabel,colorbar,axis,title,show
from matplotlib import cm

#PARAMETERS !!!!!!!!!!!!!!!!!!!!!!!!
parameters = dict(T_max = 1000, # maximal time
                  dT = 0.001, # Discretization time
                  sigma0 = 0.01,  #Initial standard variation of the population
                  x_mean0 = 0.,
                  C = 0.5,    # competition
#                 p = 0.03,      # Probability of mutation
                  b_r = 1,     # birth rate
                  d_r = 1,      # death rate
                  d_e = 2,   #exponetial power
                  beta = 0,
                  mu = 1,
                  sigma = 0.01,
                  tau = 0.3,  # transfer rate
                  X_min = -0.2, #length of the numerical interval of traits (for PDE!)
                  X_max=1.5,
                  dX = 0.01, #discretization of the space of traits
                  eps = 1
                  )
dX, T_max, dT = parameters['dX'], parameters['T_max'], parameters['dT']
X_min, X_max= parameters['X_min'], parameters['X_max']

#GRID !!!!!!!!!!!!!!!!!!

#TIME
nT = int(T_max/dT)#number of times
T = np.arange(0,T_max,dT)#space of time

nX = int((X_max-X_min)/dX) # number of traits
I=range(nX)#all the indexes
X = np.arange(X_min,X_max,dX)#space of traits

#X_large=np.arange(2*X_min,2*X_max,dX)

#INITIALIZATION !!!!!!!!!!!!!!!!!!!!!!
f0 = np.exp(-np.power(np.absolute(X-parameters['x_mean0']),2)/parameters['sigma0']*parameters['eps'])# initial density
rho0=(parameters['b_r'])/(parameters['C'])
f0=f0/(np.sum(f0)*parameters['dX'])*rho0

f = np.empty([nT, nX])#densities for all times
f[0]=f0

#LOOP OVER PARAMETERS !!!!!!!!!!!!!!!!!!!
param='tau'
param_m=0.2
param_M=0.55
param_step=0.025
J=np.arange(param_m,param_M,param_step)
for j in J:
    f[0]=f0
    parameters[param]=j
    print(' !!!!!!!!!!!!!! '+str(param)+' = '+str(j), flush = True)

    #AUXILIARY FUNCTIONS

    #DEATH
    Death= parameters['d_r']*np.power(np.absolute(X),parameters['d_e'])
    # ht_kernel = np.vectorize(lambda x: np.dot(parameters['tau'],np.heaviside(x-X, 1)))(X)

    #MUTATION
    Mutation_kernel=  np.exp(-(X/(parameters['sigma']*parameters['eps']))**2*1/2)*parameters['dX']
    constant= np.sum(Mutation_kernel)
    Mutation_kernel=parameters['b_r']*Mutation_kernel/constant
    X_min_new=int(-X_min/dX)#Bounds for the new indexes that must be kept after the convolution
    X_max_new=int((-2*X_min+X_max)/dX)

    def Next_Generation_PDE(f,parameters):
        rho = np.sum(f)*parameters['dX']
        if rho==0:
            return 0
        death_term = Death + rho*parameters['C']
        birth_part = np.convolve(Mutation_kernel,f)[X_min_new:X_max_new]
        transfer_part = parameters['tau']/rho*parameters['dX']*np.fromiter((np.sum(f[:i])-np.sum(f[i:]) for i in I),float)
        #new_f = np.maximum(0,f + dT/parameters['eps']*(np.multiply(f,(-death_term+ transfer_part))+birth_part))
        new_f = f + dT/parameters['eps']*(np.multiply(f,(-death_term+ transfer_part))+birth_part)
        return new_f

    #SIMULATION !!!!!!!!!!!!!!!!!!

    for i in range(nT-1):
        f[i+1] = Next_Generation_PDE(f[i],parameters)
        if i%1000==0:
            print('T= '+str(i*parameters['dT']), flush = True)

#PLOT !!!!!!!!!!!!!!!!!

    figure = plt.figure()
    im = imshow(f[100:].transpose()[::-1],cmap=cm.coolwarm,aspect='auto',extent=(0,parameters['T_max'],X_min,X_max),vmin=0,vmax=0.8)
    colorbar(im)

    plt.xlabel('time')
    plt.ylabel('trait')

    par_str = '' # create the title, with values of parameters
    for k, v in parameters.items():
        if k == 'N0' or k == 'd_r' or k=='X_min':
            smth = ",\n"
        else:
            smth = ", "
        par_str += k + "=" + str(v) + smth
    par_str=par_str+' rho0='+str(rho0)
    plt.title(par_str)

    plt.tight_layout()
    #levellines=np.array([-0.01])
    #cset = contour(u.transpose(),levellines,linewidths=2,cmap=cm.Set2)
    #clabel(cset,inline=True,fmt='%1.1f',fontsize=10)

    current_time = datetime.now().time()
    figure.savefig(str("horizontal_gene_transfer/Mesocentre/Figures/PDE/plot_" + str(current_time)[0:8]+".pdf")) # Possibly different delimeter on Linux and Windows!

    gc.collect()

ext: py

reason: import numpy


text:
'''
Created on May 28, 2010

@author: Shreyas Joshi
@contact: shreyasj@gatech.edu
'''

import DataAccess as da
import tables as pt
import math
import numpy as np

class Optimizer(object):
    '''
    @summary: The optimizer class is supposed to get the alpha values and data on the current portfolio to make a decision on what trades to make.

    '''

# one day in unix time

    def __init__(self, listOfStocks):
        '''
        Constructor
        '''
#        self.alphah5f= pt.openFile("randomAlpha.h5", mode = "a") # if mode ='w' is used here then the file gets overwritten!
        self.listOfLoosingStocks=list()
        self.noOfDaysStockHasBeenLoosingValue=list()

        self.listOfStocks= list(listOfStocks)
        self.DAY = 86400
        staticDataItemsList= list()
        staticDataItemsList.append("blah")
        dataItemsList= list()
        dataItemsList.append("alphaValue")
        self.minCom= 5.00
        self.ComPerShare = 0.01
#        self.alphaData= da.DataAccess(False,"randomAlpha.h5","/alphaData", "alphaData", True, listOfStocks, None, dataItemsList)
    #def __init__ ends

    def strategyOne (self, portfolio,positions,timestamp,stockInfo, dataAccess):
     output=[]
     adjOpenData= dataAccess.getMatrix (self.listOfStocks, "adj_open", timestamp- 2*self.DAY, timestamp- self.DAY)

     if adjOpenData is not None:
         #choose the biggest loser
         ctr=0
         currentBiggestLoss= - float("infinity")
         currentBiggestLoserIndex=-1
         while (ctr<len(self.listOfStocks)):
#             print "adjOpenData shape: " + str(adjOpenData.shape) + ", ctr: " + str(ctr) + ", len(self.listOfStocks): " + str(len(self.listOfStocks))
             if (adjOpenData[0][ctr] > adjOpenData[1][ctr]):
                 #Which means the stock lost value
                 if ((adjOpenData[0][ctr] - adjOpenData[1][ctr])/adjOpenData[0][ctr] > currentBiggestLoss): #biggest % loss
                     currentBiggestLoss= (adjOpenData[0][ctr] - adjOpenData[1][ctr])
                     currentBiggestLoserIndex= ctr
             ctr+=1
             #While loop done
         #Now  we have the stock which lost the most value. We buy it and also put in a sell order for 2 days later

         if (currentBiggestLoserIndex != -1):
             order= stockInfo.OutputOrder()
             order.symbol= self.listOfStocks[currentBiggestLoserIndex]
             order.volume= 10
             order.task= 'buy'
             order.orderType = 'moc'
             order.duration = self.DAY

             newOrder = order.getOutput()

             if newOrder != None:
                output.append(newOrder)
             else:
                 print "ERROR! ERROR! ERROR!"
             #if (currentBiggestLoserIndex != -1): ends
         #if adjOpenData is not None: ends
     else:
         print "adjOpenData is None!"
         #else ends
     #Now to decide which stocks to sell
     currentPositions= positions.getPositions()
     for stock in portfolio.currStocks:
         for pos in currentPositions:
          if (str(pos['symbol'])== str(stock)):
           if ((pos['timestamp'] )< timestamp - 2*self.DAY):
               temp=  dataAccess.getStockDataItem(pos['symbol'], 'adj_close', timestamp)
               if not (np.isnan(temp)):
                if ((pos['purchase_price'] + (pos['shares']*self.ComPerShare))< temp):
                 order= stockInfo.OutputOrder()
                 order.symbol= stock
                 order.volume= pos['shares']
                 order.task= 'sell'
                 order.orderType = 'moc'
                 order.closeType = 'fifo'
                 order.duration = self.DAY
                 newOrder = order.getOutput()

                 if newOrder != None:
                   output.append(newOrder)
                 else:
                   print "ERROR! ERROR! ERROR!"
        #for pos in currentPositions: ends
     return output

    def strategyTwo (self, portfolio,positions,timestamp,stockInfo, dataAccess):
        #Here we track all the stocks that continuously loose value- then buy them when they stop loosing value. Then hold them until they keep
        #gaining value. Then sell them

        output=[]

        #adjOpenData= dataAccess.getMatrix (self.listOfStocks, "adj_open", timestamp- 2*self.DAY, timestamp- self.DAY)

        adjOpenData= dataAccess.getMatrixFromTS (self.listOfStocks, "adj_open", timestamp, -1)

#        print "list of loosing stocks: "+ str (self.listOfLoosingStocks)
#        print "current positions: " + str(positions.getPositions())
#        print "no of days: " + str (self.noOfDaysStockHasBeenLoosingValue)

        if (adjOpenData is not None):
            ctr=0
            while (ctr< len(self.listOfStocks)):
                if (adjOpenData[0][ctr] > adjOpenData[1][ctr]):

                    try:
                        index2= self.listOfLoosingStocks.index(self.listOfStocks[ctr])
                        self.noOfDaysStockHasBeenLoosingValue[index2]+=1
                    except:
                        #stock not found in the list
                        self.listOfLoosingStocks.append(self.listOfStocks[ctr])
                        self.noOfDaysStockHasBeenLoosingValue.append(1)

                    currentPositions= positions.getPositions()
                    for pos in currentPositions:
                        try:
                            index2= self.listOfLoosingStocks.index(pos['symbol']) #if it isn't in this list then we don't have to sell it
                            if (self.noOfDaysStockHasBeenLoosingValue[index2] > 2):
                                # we have this stock and it lost value twice
                                # Ergo- we sell
                                #sell
                                #if ((pos['purchase_price'] + (pos['shares']*self.ComPerShare))< temp): #rig it to make money
                                  print str(pos['symbol'])+" finally lost value for "+ str(self.noOfDaysStockHasBeenLoosingValue[index2])+" days. Selling it"
                                  order= stockInfo.OutputOrder()
                                  order.symbol= pos['symbol']
                                  order.volume= pos['shares']
                                  order.task= 'sell'
                                  order.orderType = 'moc'
                                  order.closeType = 'fifo'
                                  order.duration = self.DAY

                                  newOrder = order.getOutput()
                                  if newOrder != None:
                                    output.append(newOrder)
                                  else:
                                    print "ERROR! ERROR! ERROR!"
                        except ValueError:
                            pass #index not found
                    #for pos in currentPositions

                else:
                    #this stock did not loose value

                    #Check if had been loosing value
                    #print str(self.listOfStocks[ctr])+ " gained value"

                    try:
                      index1= self.listOfLoosingStocks.index(self.listOfStocks[ctr])
                      if (self.noOfDaysStockHasBeenLoosingValue[index1]>3):

                            #print "This stock has been loosing value for atleast 3 days"

                            order= stockInfo.OutputOrder()
                            order.symbol= self.listOfStocks[ctr]
                            order.volume= min(10* self.noOfDaysStockHasBeenLoosingValue[index1], 100)
                            order.task= 'buy'
                            order.orderType = 'moc'
                            order.closeType = 'fifo'
                            order.duration = self.DAY
                            newOrder = order.getOutput()
                            if newOrder != None:
                                 output.append(newOrder)
                            else:
                                 print "ERROR! ERROR! ERROR!"

                          #The stock was loosing value for <=3 days  but now gained value- so off with the head
                      self.listOfLoosingStocks.pop(index1)
                      self.noOfDaysStockHasBeenLoosingValue.pop(index1)

                    except ValueError:
                        pass

#                    try:
#                        index1= self.listOfLoosingStocks.index(self.listOfStocks[ctr])
#
#                        print str(self.listOfStocks[ctr])+" lost value for "+ str(self.noOfDaysStockHasBeenLoosingValue[index]+ "..and then gained..")
#                        #Stock found
#                        #if it had lost value for more than 2 days then buy!
#                        if (self.noOfDaysStockHasBeenLoosingValue[index1]>3):
#                            #buy
#                            order= stockInfo.OutputOrder()
#                            order.symbol= self.listOfStocks[ctr]
#                            order.volume= max(10* self.noOfDaysStockHasBeenLoosingValue[index1], 100)
#                            order.task= 'buy'
#                            order.orderType = 'moc'
#                            order.closeType = 'fifo'
#                            order.duration = self.DAY
#                            newOrder = order.getOutput()
#                            if newOrder != None:
#                                 output.append(newOrder)
#                            else:
#                                 print "ERROR! ERROR! ERROR!"
#                        else:
#                            #it was loosing value- but for less than 2 days. So we just remove this entry...
#                            self.listOfLoosingStocks.pop(index1)
#                            self.noOfDaysStockHasBeenLoosingValue.pop(index1)
#                    except:
#                        #Not found- this stock had not lost value
#                        print "could not find index! Possibly a bug"
                ctr+=1
                #while loop ends..hopefully!

        return output
        #strategyTwo ends

    def execute (self, portfolio,positions,timestamp,stockInfo, dataAccess):
     '''
     @param portfolio: The portfolio object that has symbol and value of currently held stocks.
     @param positions: Detailed info about current stock holdings.
     @param timestamp: Current simulator time stamp
     @param stockInfo: Not used anymore for dataAccess.
     @param dataAccess: a dataAccess object that will henceforth be used to access all data
     '''

     output=[]
     #output =  self.strategyOne(portfolio, positions, timestamp, stockInfo, dataAccess)
     output =  self.strategyTwo(portfolio, positions, timestamp, stockInfo, dataAccess)
        #for pos in currentPositions: ends
     #print "The outout is: " + str(output)

     return output

#
#     #Right now this is stratDemo firstStrategy
#     output = []
#    #This first for loop goes over all of the stock data to determine which stocks to buy
#     for stock in dataAccess.getListOfStocks(): #stockInfo.getStocks(startTime = timestamp - self.DAY,endTime = timestamp):
#        # if close is higher than open and close is closer to high than open is to low, buy
#
##        print "In Optimizer"
##        print "     timestamp asked for is: " + str(timestamp - self.DAY, timestamp)
##        print "self.DAY: " + str(self.DAY)
#        adj_open= dataAccess.getStockDataList(stock, 'adj_open', timestamp - self.DAY, timestamp)
#        adj_close= dataAccess.getStockDataList(stock, 'adj_close', timestamp - self.DAY, timestamp)
#        adj_high= dataAccess.getStockDataList(stock, 'adj_high', timestamp - self.DAY, timestamp)
##        alphaValue= self.alphaData.getStockDataList (stock, 'alphaValue', timestamp - self.DAY, timestamp)
#
#        if (adj_open.size > 0):
#            #if alphaValue <= 0.5 and adj_open < adj_close and (adj_high - adj_close) > (adj_open - adj_close): #highly possible bug here?
#         if adj_open < adj_close and (adj_high - adj_close) > (adj_open - adj_close): #highly possible bug here?
#            order = stockInfo.OutputOrder()
#            order.symbol = stock #stock['symbol']
#            order.volume = 20
#            order.task = 'buy'
#            order.orderType = 'moc'
#            order.duration = self.DAY * 2
#            newOrder = order.getOutput()
#            if newOrder != None:
#                output.append(newOrder)
#
#    #This for loop goes over all of our current stocks to determine which stocks to sell
#     for stock in portfolio.currStocks:
#        openPrice = list(dataAccess.getStockDataList(stock, 'adj_open',timestamp - self.DAY, timestamp))#dataAccess.getDataList(stock, 'adj_open',timestamp - self.DAY, timestamp)#stockInfo.getPrices(timestamp - self.DAY, timestamp,stock,'adj_open')
#        closePrice = list(dataAccess.getStockDataList(stock, 'adj_close',timestamp - self.DAY, timestamp))#dataAccess.getDataList(stock, 'adj_close',timestamp - self.DAY, timestamp) #stockInfo.getPrices(timestamp - self.DAY, timestamp,stock,'adj_close')
#        highPrice =  list(dataAccess.getStockDataList(stock, 'adj_high',timestamp - self.DAY, timestamp))#dataAccess.getDataList(stock, 'adj_high',timestamp - self.DAY, timestamp) #stockInfo.getPrices(timestamp - self.DAY, timestamp,stock,'adj_high')
#        lowPrice = list(dataAccess.getStockDataList(stock, 'adj_low',timestamp - self.DAY, timestamp))#dataAccess.getDataList(stock, 'adj_low',timestamp - self.DAY, timestamp) #stockInfo.getPrices(timestamp - self.DAY, timestamp,stock,'adj_low')
#        if(len(openPrice) != 0 and len(closePrice) != 0 and len(highPrice) != 0 and len(lowPrice) != 0):
#            # if closeprice is closer to low than openprice is to high, sell
#            if (closePrice[0]-lowPrice[0]) > (highPrice[0]-openPrice[0]):
#                order = stockInfo.OutputOrder()
#                order.symbol = stock
#                order.volume = portfolio.currStocks[stock]/2+1
#                order.task = 'sell'
#                order.orderType = 'moo'
#                order.closeType = 'fifo'
#                order.duration = self.DAY * 2
#                newOrder = order.getOutput()
#                if newOrder != None:
#                    output.append(newOrder)
#    # return the sell orders and buy orders to the simulator to execute
#     return output

     #return orders
     #def execute ends

max_stars_repo_name: elxavicio/QSTK

max_issues_repo_name: jenniyanjie/QuantSoftwareToolkit

max_forks_repo_name: jenniyanjie/QuantSoftwareToolkit


text:
"""
DOCUMENTATION

This script aims to tune your model with different settings which includes randomized search, grid search or optuna

"""

from cv2 import grabCut
import numpy as np
from sklearn.ensemble import RandomForestClassifier, RandomForestRegressor
from sklearn.model_selection import GridSearchCV, RandomizedSearchCV
from xgboost import XGBRegressor

class autoTuning:

    def __init__(self, model) -> None:
        self.model = model

    def autoTuningRFr(self, train, labels):
        # n_jobs = -1 uses all the cores available
        rf = RandomForestRegressor(random_state= 42, n_jobs= -1 , verbose = 1)

        # number of trees in random forests
        n_estimators = [int(x) for x in np.linspace(start = 200 , stop = 2000, num = 10)]

        # number of features to consider at every split
        max_features = ['auto', 'sqrt']

        # maximum number of levels in tree
        max_depth = [int(x) for x in np.linspace(10, 110, num = 11)]
        max_depth.append(None)

        # min number of  samples required to split a node
        min_samples_split = [2, 5, 10]

        # min number of samples required at each leaf node
        min_samples_leaf = [1, 2, 4]

        # method of selecting samples for training each tree
        bootstrap = [True, False]

        # create a random grid

        random_grid = {

            'n_estimators' : n_estimators,
            'max_features' : max_features,
            'max_depth' : max_depth,
            'min_samples_split' : min_samples_split,
            'min_samples_leaf' : min_samples_leaf,
            'bootstrap' : bootstrap
        }

        rf_random = RandomizedSearchCV(estimator= rf, param_distributions= random_grid,
                                       n_iter= 100, cv = 3, verbose = 2, random_state=42,
                                       n_jobs= -1)

        rf_random.fit(train, labels)

        return rf_random

    def autoTuningXGBr(self, train, labels):
        """
        Fine tuning for XGBRegressor, you may increase your CV, if you have smaller
        sized dataset
        """

        model = XGBRegressor()

        parameters = {
            'nthread' : [4],
            'objective' : ['reg:linear'],
            'learning_rate' : [0.3, 0.05, .07],
            'max_depth' : [5,6,7],
            'min_child_weight' : [4],
            'silent' : [1],
            'subsample' : [0.7],
            'colsample_bytree' : [0.7],
            'n_estimators' : [350, 500, 700]
        }

        grid = GridSearchCV(model, parameters, cv = 2, n_jobs=-1, verbose = True)
        grid.fit(train, labels)
        return grid


text:
"""Functions for converting between astronomical reference frames, etc."""

import numpy as np
from casacore.measures import measures
from casacore.quanta import quantity
import math

def convertBasis(me, rbasis, from_refFrame, to_refFrame):
    basis = np.zeros((3,3))
    for comp in range(3):
        vr = np.squeeze(rbasis[:,comp])
        (az, el) = crt2sph(vr)
        vr_sph_me = measures().direction(from_refFrame,
                                        quantity(az, 'rad'),
                                        quantity(el, 'rad'))
        v_sph_me = me.measure(vr_sph_me, to_refFrame)
        v_me = sph2crt_me(v_sph_me)
        basis[:,comp] = v_me
    return basis

def CEL2TOPOpnts(obsTimes, stnPos, celPnt):
    #Convert python times to pyrap times
    obsTimes_lst = []
    for obsTime in obsTimes:
        obsTimes_lst.append(quantity(obsTime.isoformat()).get_value())
    obsTimes_me = quantity(obsTimes_lst,'d')
    #Convert source direction to pyrap
    celPntTheta, celPntPhi, celPntRefFrame = celPnt
    celPnt = celPntRefFrame, str(celPntPhi), str(math.pi/2-celPntTheta)
    celPnt_me = measures().direction(celPnt[0],
                                   celPnt[1]+'rad',
                                   celPnt[2]+'rad')
    stnPos_me = measures().position('ITRF', str(stnPos[0,0])+'m',
                            str(stnPos[1,0])+'m',
                            str(stnPos[2,0])+'m')
    celPntBasis = getSph2CartTransf(sph2crt_me(celPnt_me))
    obsTimesArr = obsTimes_me.get_value()
    obsTimeUnit = obsTimes_me.get_unit()
    #CelRot=zeros((len(obsTimesArr),2,2))
    rotang = np.zeros(len(obsTimesArr))
    me = measures()
    #Set position of reference frame w.r.t. ITRF
    me.doframe(stnPos_me)

    me.doframe(me.epoch('UTC', quantity(obsTimesArr[0], obsTimeUnit)))
    CelRot0 = getRotbetweenRefFrames(celPntRefFrame,'ITRF', me)
    rotang[0] = 0.0

    for ti in range(1,len(obsTimesArr)):
        #Set current time in reference frame
        timEpoch = me.epoch('UTC',quantity(obsTimesArr[ti],obsTimeUnit))
        me.doframe(timEpoch)
        #Incomplete
        CelRot = getRotbetweenRefFrames(celPntRefFrame,'ITRF', me)
        IncRot = CelRot*CelRot0.T
        rotang[ti] = rotzMat2ang(IncRot)
    return CelRot0, rotang

def getParallacticRot(obsTimes, stnPos, srcDir, doPolPrec=True):
    #Convert python times to pyrap times
    obsTimes_lst = []
    for obsTime in obsTimes:
        obsTimes_lst.append(quantity(obsTime.isoformat()).get_value())
    obsTimes_me = quantity(obsTimes_lst,'d')
    print obsTimes_me
    #Convert source direction to pyrap
    srcTheta, srcPhi, srcRefFrame = srcDir
    srcDir = srcRefFrame, str(srcPhi), str(math.pi/2-srcTheta)
    srcDir_me = measures().direction(srcDir[0],
                                   srcDir[1]+'rad',
                                   srcDir[2]+'rad')
    stnPos_me = measures().position('ITRF',str(stnPos[0,0])+'m',
                            str(stnPos[1,0])+'m',
                            str(stnPos[2,0])+'m')

    obsTimesArr = obsTimes_me.get_value()
    obsTimeUnit = obsTimes_me.get_unit()
    paraMat = np.zeros((len(obsTimesArr), 2, 2))

    me = measures()
    #Set position of reference frame w.r.t. ITRF
    me.doframe(stnPos_me)

    if doPolPrec:
        #Get sky precession rotation matrix
        #(Assuming no change over data interval)
        me.doframe(me.epoch('UTC', quantity(obsTimesArr[0], obsTimeUnit)))
        precMat = getSkyPrecessionMat(me,srcDir_me)
    for ti in range(len(obsTimesArr)):
        #Set current time in reference frame
        timEpoch = me.epoch('UTC', quantity(obsTimesArr[ti], obsTimeUnit))
        me.doframe(timEpoch)

        #Compute polariz comps in spherical sys to cartesian Station coord sys
        #paraMtc=computeParaMat_tc('J2000', 'ITRF', srcDir_me, me)

        #Alternatively:
        paraMme = computeParaMat_me('J2000', 'AZEL', srcDir_me, me)

        paraM = paraMme

        if doPolPrec:
            #With precession:
            paraM=paraM * precMat
            #else:
            #Do not apply precession rotation of polarimetric frame.
            #This is then the apparent polarization frame.

        paraMat[ti,:,:] = paraM
    return paraMat

def getSkyPrecessionMat(me, srcDirection):
    """Compute precession matrix. This is the 2D rotation along srcDirection
    between J2000 and current epoch.

    At J2000 epoch:
    compute 2 cartesian vectors orthogonal to srcDirection: alpha & delta.
    alpha is the orthogonal direction on equator."""
    alphaJ2000ra = srcDirection['m0']['value']+math.pi/2
    alphaJ2000dec = 0.0
    alphaJ2000vec = sph2crt(alphaJ2000ra,alphaJ2000dec)
    #delta is the orthogonal direction along the meridian.
    deltaJ2000ra = srcDirection['m0']['value']
    deltaJ2000dec = srcDirection['m1']['value']+math.pi/2
    deltaJ2000vec = sph2crt(deltaJ2000ra,deltaJ2000dec)

    alpha = me.direction('J2000',
               str(alphaJ2000ra)+'rad',str(alphaJ2000dec)+'rad')
    delta = me.direction('J2000',
               str(deltaJ2000ra)+'rad',str(deltaJ2000dec)+'rad')
    #Convert alpha & delta to directions in the current epoch
    alphaTru = me.measure(alpha,'JTRUE') #'JTRUE' isn't stable
    deltaTru = me.measure(delta,'JTRUE')
    raA = alphaTru['m0']['value']
    decA = alphaTru['m1']['value']
    alphaTruvec = sph2crt(raA,decA)
    raD = deltaTru['m0']['value']
    decD = deltaTru['m1']['value']
    deltaTruvec = sph2crt(raD, decD)

    cosPrecAng = ((alphaTruvec*alphaJ2000vec.T)[0,0]
                +(deltaTruvec*deltaJ2000vec.T)[0,0])/2.0
    sinPrecAng = ((alphaTruvec*deltaJ2000vec.T)[0,0]
                -(deltaTruvec*alphaJ2000vec.T)[0,0])/2.0
    #Precession of polarization basis is
    #precMat=np.matrix([ [(alphaTruvec*alphaJ2000vec.T)[0,0],
    #                     (alphaTruvec*deltaJ2000vec.T)[0,0]],
    #                    [(deltaTruvec*alphaJ2000vec.T)[0,0],
    #                     (deltaTruvec*deltaJ2000vec.T)[0,0]] ])
    precMat = np.matrix([[ cosPrecAng,sinPrecAng],
                         [-sinPrecAng,cosPrecAng]])
    #print precMat
    return precMat

def getRotbetweenRefFrames(rfFrom,rfTo, me):
    x_fr = me.direction(rfFrom,  '0deg',  '0deg')
    y_fr = me.direction(rfFrom, '90deg',  '0deg')
    z_fr = me.direction(rfFrom, '90deg', '90deg')
    x_to = np.asmatrix(sph2crt_me(me.direction(rfTo,    '0deg',  '0deg')))
    y_to = np.asmatrix(sph2crt_me(me.direction(rfTo,   '90deg',  '0deg')))
    z_to = np.asmatrix(sph2crt_me(me.direction(rfTo,   '90deg', '90deg')))
    x_fr_to = np.asmatrix(sph2crt_me(me.measure(x_fr, rfTo)))
    y_fr_to = np.asmatrix(sph2crt_me(me.measure(y_fr, rfTo)))
    z_fr_to = np.asmatrix(sph2crt_me(me.measure(z_fr, rfTo)))
    xyz_fr_to = np.bmat([[x_fr_to],[y_fr_to],[z_fr_to]])
    xyz_to = np.bmat([[x_to],[y_to],[z_to]])
    frameRot = xyz_fr_to*xyz_to.T
    return frameRot

def computeParaMat_me(rfFrom, rfTo, aDir_me, me):
    """Compute parallactic rotation matrix. me contains epoch"""
    NEvFrom = computeSphBasis(rfFrom, rfTo, 0, aDir_me, me)
    NEvTo = computeSphBasis(rfFrom, rfTo, 1, aDir_me, me)
    paraMat = NEvTo * NEvFrom.T
    return paraMat

def computeSphBasis(rfFrom, rfTo, order, aDir_me, me):
    RbaseSph = me.measure(aDir_me, rfTo)
    Rbase = sph2crt_me(RbaseSph)
    if order == 0:
        NbaseFrom = me.direction(rfFrom,
                         str(aDir_me['m0']['value'])+'rad',
                         str(aDir_me['m1']['value']+math.pi/2.0)+'rad'
                        )
        NbaseSph = me.measure(NbaseFrom, rfTo)
    else:
        NbaseSph = me.direction(rfTo,
                         str(RbaseSph['m0']['value'])+'rad',
                         str(RbaseSph['m1']['value']+math.pi/2.0)+'rad'
                        )
    Nbase = sph2crt_me(NbaseSph)
    Ebase = np.cross(Nbase, Rbase)
    NEbasis = np.bmat([[Nbase],[Ebase]])
    return NEbasis

def computeParaMat_tc(rfFrom, rfTo, aDir_me, me):
    framRot=getRotbetweenRefFrames(rfFrom, rfTo, me)
    srcFromcrt=sph2crt_me(me.measure(aDir_me,rfFrom))
    pol2cartFrom=computeSph2CrtMat(srcFromcrt)
    srcTocrt=sph2crt_me(me.measure(aDir_me,rfTo))
    pol2cartTo=computeSph2CrtMat(srcTocrt)
    paraMat=pol2cartFrom.T*framRot*pol2cartTo
    return paraMat

def getSph2CartTransf(r):
    """Compute the transformation matrix from a spherical basis to a Cartesian
    basis at the field point given by the input 'r'.
    The output 'transf_sph2cart' is defined such that:

    [[v_x], [v_y], [v_z]]=transf_sph2cart*matrix([[v_r], [v_phi], [v_theta]])."""
    r = np.matrix(r)
    r = r.squeeze()
    ru = r/np.sqrt((r*r.T)[0,0])
    xu = ru[0,0]
    yu = ru[0,1]
    zu = ru[0,2]
    tang2cart = 1.0/np.sqrt(xu*xu+yu*yu)* np.matrix([
                          [ yu,          xu*zu],
                          [-xu,          yu*zu],
                          [ 0., -(xu*xu+yu*yu)]])
    transf_sph2cart = np.bmat([[ru.T, tang2cart]])
    return transf_sph2cart

def getSph2CartTransfArr(r):
    """(Array version) Compute the transformation matrix from a spherical basis to a Cartesian
    basis at the field point given by the input 'r'.
    The output 'transf_sph2cart' is defined such that:

    [[v_x], [v_y], [v_z]]=transf_sph2cart*matrix([[v_r], [v_phi], [v_theta]])."""
    r = np.array(r)
    r = r.squeeze()
    #ru = r/np.sqrt(np.tensordot(r,r, axes=([0],[0])))
    ru = r/np.sqrt(r[0,...]**2+r[1,...]**2+r[2,...]**2)
    xu = ru[0,...]
    yu = ru[1,...]
    zu = ru[2,...]
    nrf = 1.0/np.sqrt(xu*xu+yu*yu)
    transf_sph2cart = np.array([
                          [xu,             yu*nrf,          xu*zu*nrf],
                          [yu,            -xu*nrf,          yu*zu*nrf],
                          [zu, np.zeros(xu.shape), -(xu*xu+yu*yu)*nrf]])
    return transf_sph2cart

def computeSph2CrtMat(lmnMatrix):
    lmn = lmnMatrix.squeeze()
    l = lmn[0,0]
    m = lmn[0,1]
    n = lmn[0,2]
    polz2cart = 1.0/np.sqrt(l*l+m*m)* np.matrix([
                          [  m,        l*n],
                          [ -l,        m*n],
                          [ .0, -(l*l+m*m)]])
    return polz2cart

def rotzMat2ang(rotMat):
    ang = np.arctan2(rotMat[0,1], rotMat[0,0])
    return ang

def crt2sph(dir_crt):
    x = np.squeeze(dir_crt[0])
    y = np.squeeze(dir_crt[1])
    z = np.squeeze(dir_crt[2])
    azi = np.arctan2(y, x)
    ele = np.arcsin(z)
    #Radial component (not needed for unit directions)
    #r=np.sqrt(x**2+y**2+z**2)
    #ele=np.arcsin(z/r)
    return azi, ele

def sph2crt_me(sphme):
    return sph2crt(sphme['m0']['value'], sphme['m1']['value'])

def sph2crt(azi, ele):
    #Spherical polar angles in azimuth and elevation to cartesian conversion
    x = np.cos(ele)*np.cos(azi)
    y = np.cos(ele)*np.sin(azi)
    z = np.sin(ele)
    return(np.array([x,y,z]))
    #return(np.matrix([x,y,z]))

def IAU_pol_basis(src_az, src_el):
    """Compute the (x_hat, y_hat, z_hat) basis in IAU polarization system for a direction
    given by (azimuth, elevation) tuple typically Ra, Dec."""
    #Here C09 refers to Carozzi2009 vCZ paper: it has (r_hat, phi_hat, theta_hat).
    #While IAU has +x pointing to -theta_hat, +y along +phi_hat, +z along -r_hat.
    #HOWEVER, implementation here keeps r_hat in first columns, so:
    IAUtoC09 = np.array([[ 1.,  0.,  0.],
                         [ 0.,  0., -1.],
                         [ 0., +1.,  0.]])
    #IAUtoC09 = np.identity(3)
    basis_C09 = np.array(getSph2CartTransf(sph2crt(src_az, src_el)))
    basis_IAU = np.matmul(basis_C09, IAUtoC09)
    return basis_IAU

def pyTimes2meTimes(pyTimes):
    obsTimes_lst = []
    for obsTime in pyTimes:
        obsTimes_lst.append(quantity(obsTime.isoformat()).get_value())
    obsTimes_me = quantity(obsTimes_lst,'d')
    obsTimesArr = obsTimes_me.get_value()
    obsTimeUnit = obsTimes_me.get_unit()
    return obsTimesArr, obsTimeUnit

def setEpoch(obsTimesArr,obsTimeUnit):
    stnPos_me = measures().position('ITRF', '0m', '0m', '0m')
    me = measures()
    #Set position of reference frame w.r.t. ITRF
    me.doframe(stnPos_me)
    timEpoch = me.epoch('UTC', quantity(obsTimesArr, obsTimeUnit))
    me.doframe(timEpoch)
    return me

def printJones(Jn):
    for ti in range(0,Jn.shape[0]):
        print Jn[ti,:,:]

def shiftmat2back(arr):
    arr=np.rollaxis(arr,0,arr.ndim)
    arr=np.rollaxis(arr,0,arr.ndim)
    return arr

max_stars_repo_name: creaneroDIAS/dreamBeam

max_issues_repo_name: creaneroDIAS/dreamBeam

max_forks_repo_name: creaneroDIAS/dreamBeam


text:
import pandas as pd
import numpy as np
import sys
from keras.callbacks import ModelCheckpoint
from keras.callbacks import CSVLogger
from keras.callbacks import ReduceLROnPlateau
from keras.models import Sequential
from keras.layers import Dense, Dropout
from keras.utils import np_utils
from sklearn.preprocessing import LabelEncoder
import os
import time
import h5py
import numpy as np
from keras.callbacks import Callback
from sklearn.metrics import confusion_matrix, f1_score, precision_score, recall_score
from sklearn.externals import joblib

trainName = sys.argv[1]
testName = sys.argv[2]

networkTopologies = [["M-M", 15, 15]]

# Create an object called iris with the iris Data
dftrain = pd.read_csv(filepath_or_buffer=trainName, header=None, sep=',')
dftest = pd.read_csv(filepath_or_buffer=testName, header=None, sep=',')

cols = ['Proto']
for i in range(1,dftrain.shape[1]):
    cols.append('Byte' + str(i))

dftrain.columns=cols
dftrain.dropna(how="all", inplace=True)
dftrain.tail()

dftest.columns=cols
dftest.dropna(how="all", inplace=True)
dftest.tail()

Xtrain = dftrain.ix[:,1:dftrain.shape[1]].values
Ytrain = dftrain.ix[:,0].values
Xtest = dftest.ix[:,1:dftest.shape[1]].values
Ytest = dftest.ix[:,0].values

encoder = LabelEncoder()
encoder.fit(Ytrain)
encYtrain = encoder.transform(Ytrain)

encoder = LabelEncoder()
encoder.fit(Ytest)
encYtest = encoder.transform(Ytest)

joblib.dump(encoder, "NNENC.pkl")

filename = "models/ABD/NN/f1.csv"

class Metrics(Callback):
    def on_train_begin(self, logs={}):
        with open(filename, "w") as file:
            file.write("epoch,f1\n")

    def on_epoch_end(self, epoch, logs={}):
        val_predict = (np.asarray(self.model.predict(self.validation_data[0]))).round()
        val_targ = self.validation_data[1]
        f1 = f1_score(val_targ, val_predict)
        with open(filename, "a") as file:
            file.write(str(epoch) + "," + str(f1) + "\n")
        return

f1Cal = Metrics()

directory = "models/ABD/NN/"
if not os.path.exists(directory):
    os.makedirs(directory)

#Step into new dir
for i in xrange(0, 1):
    logfile = directory + "log-" + str(i) + ".csv"

    print("[*] STARTING NN")
    print("    Topology   : " + networkTopologies[0][0])

    model = Sequential()
    model.add(Dense(networkTopologies[0][1],activation='relu', kernel_initializer='uniform', input_shape=(dftrain.shape[1]-1,)))
    model.add(Dropout(0.5))

    for k in xrange(2, len(networkTopologies)):
        model.add(Dense(networkTopologies[0][k],activation='relu', kernel_initializer='uniform',))
        model.add(Dropout(0.5))

    model.add(Dense(1, activation='sigmoid', kernel_initializer='uniform',))
    model.compile(optimizer='rmsprop', loss='binary_crossentropy', metrics=['accuracy'])

    reduce_lr = ReduceLROnPlateau(monitor='val_loss', factor=0.2, patience=5, min_lr=0.001)
    log = CSVLogger(logfile, separator=',', append=True)
    checkpoint = ModelCheckpoint(directory + "/weights.{epoch:04d}-{val_acc:.04f}.hdf5", monitor='val_acc', verbose=0,
                                 save_best_only=True, save_weights_only=False, mode='auto', period=1)
    model.fit(Xtrain, encYtrain, batch_size=64, epochs=100, verbose=1, validation_data=(Xtest,encYtest),
              callbacks=[log, reduce_lr, f1Cal, checkpoint])

    # # evaluate the model
    # scores = model.evaluate(Xtest, encYtest, batch_size=128)
    # print("\n%s: %.2f%%" % (model.metrics_names[1], scores[1]*100))
    #
    # # Create actual english names for the plants for each predicted plant class
    # preds = [model.predict_classes(Xtest, batch_size=32)]
    # certainties = model.predict(Xtest, batch_size=32)
    #
    # numRecords = certainties.shape[0]
    # probRunner = 0.0
    # for i in certainties:
    #     probRunner += np.amax(i)
    #
    # probRunner /= numRecords
    # print("Avg certainty : " + str(probRunner))
    #
    # if(len(preds) > 0):
    #     preds = list(encoder.inverse_transform(preds))

# print(dftest['Proto'].head())
# df = pd.crosstab(dftest['Proto'], preds, rownames=['Actual Protocol'], colnames=['Predicted Protocol'])
# df.to_csv('BinaryConfusionMatrix.csv')


text:
import numpy as np
import tensorflow as tf
import sleap
from sleap.nn.data import offset_regression

sleap.use_cpu_only()  # hide GPUs for test

def test_make_offsets():
    points = np.array([[1.8, 2.3], [0.4, 3.1], [np.nan, np.nan]], "float32")
    xv, yv = sleap.nn.data.confidence_maps.make_grid_vectors(4, 4, output_stride=1)
    off = offset_regression.make_offsets(points, xv, yv, stride=1)

    assert off.shape == (4, 4, 3, 2)

    XX, YY = tf.meshgrid(xv, yv)
    np.testing.assert_allclose(XX + off[..., 0, 0], 1.8, atol=1e-6)
    np.testing.assert_allclose(YY + off[..., 0, 1], 2.3, atol=1e-6)
    np.testing.assert_allclose(XX + off[..., 1, 0], 0.4, atol=1e-6)
    np.testing.assert_allclose(YY + off[..., 1, 1], 3.1, atol=1e-6)
    np.testing.assert_array_equal(off[..., 2, :], 0)

def test_mask_offsets():
    points = np.array([[1.0, 1.0]], "float32")
    xv, yv = sleap.nn.data.confidence_maps.make_grid_vectors(4, 4, output_stride=1)
    off = offset_regression.make_offsets(points, xv, yv, stride=1)
    cm = sleap.nn.data.confidence_maps.make_confmaps(points, xv, yv, sigma=1)
    off_mask = offset_regression.mask_offsets(off, cm, threshold=0.2)

    np.testing.assert_array_equal(off_mask[:3, :3], off[:3, :3])
    np.testing.assert_array_equal(off_mask[3:, :], 0.0)
    np.testing.assert_array_equal(off_mask[:, 3:], 0.0)

max_stars_repo_name: hectorcarrion/sleap

max_issues_repo_name: oeway/sleap

max_forks_repo_name: oeway/sleap


text:
import os
import sys
os.chdir('../../')
#print(os.listdir('.'))
sys.path.append('.')
import pickle
import pandas as pd
import numpy as np
from ml.model import MLAF
from sklearn.ensemble import GradientBoostingRegressor
import time
from sklearn.model_selection import train_test_split

def relative_error(y_true, y_hat):
    return np.mean(np.abs((y_true-y_hat)/y_true))

def f_relative_error(y_true, dtrain):
    y_hat = dtrain.get_label()
    return np.mean(np.abs((y_true-y_hat)/y_true))

def gradient(preds, dtrain):
    #Gradient for custom error
    y = dtrain.get_label()
    return (2*(preds-y)) / np.power(y,2)

def hessian(preds, dtrain):
    y = dtrain.get_label()
    return 2/np.power(y,2)

def custom_relative_error_loss(preds, dtrain):
    grad = gradient(preds, dtrain)
    hess = hessian(preds, dtrain)
    return grad, hess

MODEL_CATALOGUE = {}

qdf = pd.read_pickle('input/instacart_queries/qdf.pkl')
targets = [name  for name in qdf.columns if 'lb' not in name and 'ub' not in name]
#Filtering out the product of joins in the aggregate functions
sum_columns = [name for name in qdf[targets].columns if 'sum' in name]
avg_columns = [name for name in qdf[targets].columns if 'avg' in name]
count_columns = [name for name in qdf[targets].columns if 'count' in name]
#Generate Dataframes per aggregate function
sum_df = qdf.iloc[qdf['sum_add_to_cart_order'].dropna(axis=0).index]
avg_df = qdf.iloc[qdf['avg_add_to_cart_order'].dropna(axis=0).index]
count_df = qdf.iloc[qdf['count'].dropna(axis=0).index]

features = [name for name in sum_df.columns if name not in ['sum_add_to_cart_order','avg_add_to_cart_order','count']]

target_sum = 'sum_add_to_cart_order'
target_avg = 'avg_add_to_cart_order'
target_count = 'count'
count_df = count_df[(count_df[target_count]!=0)] # Remove 0 because it produces an error on relative error
count_df['product_name_lb'] = count_df['product_name_lb'].replace(np.nan, 'isnone')
count_df['product_name_lb'] = count_df['product_name_lb'].astype('category')
labels =  count_df['product_name_lb'].cat.codes
categorical_attribute_catalogue = {key : value for key,value in zip(count_df['product_name_lb'].values, labels)}
count_df['product_name_lb'] = labels

with open('catalogues/labels_catalogue.pkl', 'wb') as f:
    pickle.dump(categorical_attribute_catalogue,f)

avg_df[target_avg] = avg_df[target_avg].astype(float)
models_train = [(sum_df, target_sum, 'sum_add_to_cart_order'), (avg_df, target_avg, 'avg_add_to_cart_order'), (count_df, target_count, 'count')]

del qdf
# # read in data
for df, label,af in models_train:
    df = df.iloc[:5000]
    X = df[features].replace(np.nan, -1).values
#    X = np.nan_to_num(X)
    y = df[label].values
    X_train, X_test, y_train, y_test = train_test_split(X, y, train_size=0.7, random_state=1234)

    clf = GradientBoostingRegressor()
    start = time.time()
    clf.fit(X_train, y_train)

    rel_error =relative_error(y_test, clf.predict(X_test))
    print("Relative Error for {} is {}".format(label, rel_error))
    print("Time to train for {} \t took : {}".format(label, time.time()-start))

    ml_est = MLAF(clf, rel_error, features, label)
    MODEL_CATALOGUE[af] = ml_est
    # xgb_model.save_model('/home/fotis/dev_projects/model-based-aqp/catalogues/{}.dict_model'.format(label))
with open('code/model-serving/model_catalogue.pkl', 'wb') as f:
    pickle.dump(MODEL_CATALOGUE, f)

max_stars_repo_name: Skeftical/modelbasedaqp

max_issues_repo_name: traindb-project/traindb-benchmark

max_forks_repo_name: traindb-project/traindb-benchmark


text:
import numpy as np
from fractions import Fraction as Q

eps = 10e-4

def prvect(v):
	print("( ", end="")
	for i in v:
		print(format(float(i), ".4f"), end=" ")
	print(")")

def prmatr(m):
	print("[")
	for i in m:
		prvect(i)
	print("]")

def swap_rows(A, x, y):
	if x == y: return A
	temp = A[x].copy()
	A[x] = A[y]
	A[y] = temp
	return(A)

def make_m(main, ind):
	m = np.eye(3, dtype=Q)
	temp = main[ind, ind]
	for i in range(ind, 3):
		m[i, ind] = - main[i, ind] / temp if i != ind else 1 / temp
	return m

def make_rational(a):
	temp = np.zeros(shape=(3, 4), dtype=Q)
	for i in range(0, 3):
		for j in range(0, 4):
			temp[i, j] = Q(a[i, j])
	return temp

def count_cordinality_num(m):
	# det = m[0, 0] * m[1, 1] * m[2, 2] + m[2, 0] * m[0, 1] * m[1, 2] + m[0, 2] * m[1, 0] * m[2, 1] - m[0, 2] * m[1, 1] * m[2, 0] -  m[1, 0] * m[0, 1] * m[2, 2] - m[2, 1] * m[1, 2] * m[0, 0]
	inv = np.linalg.inv(raw_main_matrix)
	return np.linalg.norm(inv, np.inf) * np.linalg.norm(raw_main_matrix, np.inf)

def is_correct():
	print("DET = {}".format(format(np.linalg.det(raw_main_matrix), ".4g")))
	print("DET > 0: ", np.linalg.det(raw_main_matrix) > 0)
	if np.linalg.det(raw_main_matrix) == 0:
		print("ERROR, INCORRECT MATRIX")
		return False
	print("\nA = A^T: ", np.array_equal(raw_main_matrix, raw_main_matrix.T))
	print("DIAGONAL PRIORITY: ", check_diagonal_priority(raw_main_matrix))
	print()
	return True

def check_diagonal_priority(m):
	for i in range(len(m)):
		sum_n = 0
		for j in range(len(m[i])):
			if i != j:
				sum_n += abs(m[i][j])
		if sum_n > m[i][i]:
			return False
	return True

def solve_Gauss(matr, const):

	P_arr = []
	M_arr = []

	for k in range(0, 3):
		abs_max_indexes = np.absolute(matr).argmax(0)
		P = np.eye(3, dtype=Q)
		P = swap_rows(P, abs_max_indexes[k], k)
		print("\nP{}: ".format(k))
		prmatr(P)
		P_arr.append(P)
		matr = P.dot(matr)
		const = P.dot(const)
		M = make_m(matr, k)
		print("\nM{}: ".format(k))
		prmatr(M)
		M_arr.append(M)
		matr = M.dot(matr)
		const = M.dot(const)
		print("\nA{}:".format(k))
		prmatr(matr)
		print("\nCONST:", end="")
		prvect(const)
		print("\n>------------------------------<\n")

	x = [0] * 3
	x[2] = const[2]
	x[1] = const[1] - matr[1,2] * x[2]
	x[0] = const[0] - matr[0,2] * x[2] - matr[0,1] * x[1]
	print("X: ", x)

	E_arr = [np.eye(3, dtype=Q)]
	for k in range(0, 3):
		E = M_arr[k].dot(P_arr[k]).dot(E_arr[k])
		print("\nE{}: ".format(k+1))
		prmatr(E)
		E_arr.append(E)

	E_f = E_arr[3]

	a33 = E_f[2, 2]
	a32 = E_f[2, 1]
	a31 = E_f[2, 0]
	a21 = E_f[1, 0] - matr[1, 2] * a31
	a22 = E_f[1, 1] - matr[1, 2] * a32
	a23 = E_f[1, 2] - matr[1, 2] * a33
	a11 = E_f[0, 0] - matr[0, 1] * a21 - matr[0, 2] * a31
	a12 = E_f[0, 1] - matr[0, 1] * a22 - matr[0, 2] * a32
	a13 = E_f[0, 2] - matr[0, 1] * a23 - matr[0, 2] * a33

	final = np.array([
		[a11, a12, a13],
		[a21, a22, a23],
		[a31, a32, a33]
	])
	print("\nFINAL EQUALITY: ")
	prmatr(matr)
	print("*")
	prmatr(final)
	print("=")
	prmatr(E_f)
	print("\nEQUALS?: ")
	print(np.array_equal(matr.dot(final), E_f))

def check_stop(x1, x2):
	if abs(x1[0] - x2[0]) < eps or abs(x1[1] - x2[1]) < eps or abs(x1[2] - x2[2]) < eps:
		return True
	return False

def solve_Jacobi(matr, const):
	const[0] = const[0] / matr[0, 0]
	matr[0] /= matr[0, 0]
	const[1] /= matr[1, 1]
	matr[1] /= matr[1, 1]
	const[2] /= matr[2, 2]
	matr[2] /= matr[2, 2]
	X = np.array([Q(0)] * 3)
	X_prev = const.copy()
	X[0] = - matr[0, 1] * X_prev[1] - matr[0, 2] * X_prev[2] + const[0]
	X[1] = - matr[1, 0] * X_prev[0] - matr[1, 2] * X_prev[2] + const[1]
	X[2] = - matr[2, 0] * X_prev[0] - matr[2, 1] * X_prev[1] + const[2]
	i = 1
	print("X0 = ", end="")
	prvect(X_prev)
	print("X1 = ", end="")
	prvect(X)
	while np.linalg.norm(X_prev - X, np.inf) > eps:
		X_prev = X.copy()
		i+=1
		X[0] = - matr[0, 1] * X_prev[1] - matr[0, 2] * X_prev[2] + const[0]
		X[1] = - matr[1, 0] * X_prev[0] - matr[1, 2] * X_prev[2] + const[1]
		X[2] = - matr[2, 0] * X_prev[0] - matr[2, 1] * X_prev[1] + const[2]
		print("X{} = ".format(i), end="")
		prvect(X)

	print("RESULT X: ", end="")
	prvect(X)

if __name__ == "__main__":
	raw = np.loadtxt('matrix.txt')
	raw_main_matrix = raw[:,:-1]
	raw = make_rational(raw)
	main_matrix = raw[:,:-1]
	constants = raw[:,-1]
	if is_correct():
		print("Cond number:")
		print(count_cordinality_num(main_matrix))
		print("\nGauss method:")
		solve_Gauss(main_matrix.copy(), constants.copy())
		print("\n>------------------------------<\n\n>------------------------------<\n\nJacobi method:\n")
		solve_Jacobi(main_matrix.copy(), constants.copy())