text:
from __future__ import division
import math
import numpy as np
import random
from graph import Graph, Edge
from search_classes import SearchNode, Path
from utils import *
from shapely.geometry import Point, LineString

class PRMPathPlanner():

    """Plans path using PRM algorithm on a 2D environment.

    """
    def initialise(self, environment, bounds, start_pose, goal_region, radius, resolution, isLazy):
        """
        Initialises the main parameters for the PRM based path planner.

        Args:
            environment (a yaml environment): Environment where the planner will run.
            bounds ( (int int int int) ): min x, min y, max x, and max y coordinates of the bounds of the world.
            radius (float): Radius of the object
            resolution (int): Number of segments used to approximate a circle around a point
            isLazy (bool): If true, graph is created and collisions are only detected after running A* to find the optimum path instead of during the creation of the graph.
            goal_region: (Polygon): A polygon representing the region we want our object to reach.
            start_pose: ( ( float float) ): Starting x and y coordinates of the object in question.

        Returns:
            None
        """
        self.env = environment
        self.bounds = bounds
        self.radius = radius
        self.resolution = resolution
        self.isLazy = isLazy
        self.goal_region = goal_region
        self.start_pose = start_pose

    def path (self,environment, bounds, q_init_point, q_goal_region, radius, resolution, isLazy):
        """
        Creates and returns a path along with the set of vertices and edges that make up the graph used in the PRM algorithm.

        Args:
            environment (a yaml environment): Environment where the planner will run.
            bounds ( (int int int int) ): min x, min y, max x, and max y coordinates of the bounds of the world.
            radius (float): Radius of the object
            resolution (int): Number of segments used to approximate a circle around a point
            isLazy (bool): If true, graph is created and collisions are only detected after running A* to find the optimum path instead of during the creation of the graph.
            goal_region: (Polygon): A polygon representing the region we want our object to reach.
            start_pose: ( ( float float) ): Starting x and y coordinates of the object in question.

        Returns:
            path (list<(int,int)>): A list of tuples/coordinates representing the nodes in a path from start to the goal region
            self.V (set<(int,int)>): A set of Vertices (coordinates) of nodes in the tree
            self.E (set<(int,int),(int,int)>): A set of Edges connecting one node to another node in the tree
        """
        self.initialise(environment, bounds, q_init_point, q_goal_region, radius, resolution, isLazy)

        x0, y0 = q_init_point
        x1, y1 = q_goal_region.centroid.coords[0]
        q_init = (x0, y0)
        q_goal = (x1, y1)
        V = set()
        E = set()

        if q_init == q_goal:
            path = [q_init, q_goal]
            V.union([q_init, q_goall])
            E.union([(q_init, q_goal)])
        elif self.isEdgeCollisionFree(q_init, q_goal, self.radius):
            path = [q_init, q_goal]
            V.union([q_init, q_goal])
            E.union([(q_init, q_goal)])
        else:
            path, V, E = self.PRMSearch(self.bounds, q_init, self.radius, q_goal, isLazy)

        return path, V, E

    def PRMSearch(self, bounds, q_init, radius, q_goal, isLazy):
        '''print "starting roadmap construction"'''
        # number of nodes to put in the roadmap
        n = min(max(len(self.env.obstacles) *2, 500), 2000)
        # number of closest neighbours to examine for each configuration
        k = min(int(n/50),15)

        # construct the probabilistic roadmap
        V,E = self.roadmapConstruction(bounds, radius, k, n, isLazy)

        # k closest neighbours of start configuration q_init, from V, sorted by distance
        Nqinit = self.find_k_closest_neighbours(V, q_init, k)
        # k closest neighbours of goal configuration q_goal, from V, sorted by distance
        Nqgoal = self.find_k_closest_neighbours(V, q_goal, k)
        # V is the union of itself and Nqinit and Nqgoal
        V = V.union(Nqinit, Nqgoal)

        # iterate through Nqinit and try to find the closest neighbour to q_init that has a collision free edge between them
        for qprime in Nqinit:
            if self.isEdgeCollisionFree(q_init, qprime, radius):
                E.add((q_init, qprime))
                break

        # iterate through Nqgoal and try to find the closest neighbour to q_goal that has a collision free edge between them
        for qprime in Nqgoal:
            if self.isEdgeCollisionFree(q_goal, qprime, radius):
                E.add((q_goal, qprime))
                break
        '''print "starting search"'''
        # find, using A* Search, the shortest path between our intial point and goal point, given the constructed graph.
        path = None
        if isLazy:
            ok_edges = set()
            attempt = 1

            while True:

                path_found = True
                path = self.findShortestPath(V, E, q_init, q_goal)

                if path == None:
                    if attempt == 5:
                        break
                    '''print attempt'''
                    n += 200
                    V, E = self.roadmapReconstruction(bounds, radius, k, n, isLazy, V, E, q_init, q_goal)
                    attempt += 1
                    continue

                for i in xrange(len(path)-1):
                    edge = (path[i], path[i+1])
                    reversed_edge = (path[i+1], path[i])

                    if edge not in ok_edges:
                        if not self.isEdgeCollisionFree(path[i], path[i+1], radius):
                            if edge in E:
                                E.remove(edge)
                            if reversed_edge in E:
                                E.remove(reversed_edge)
                            path_found = False
                        else:
                            ok_edges.add(edge)
                            ok_edges.add(reversed_edge)
                if path_found:
                    break

        else:
            for attempt in xrange(1,6):
                '''print attempt'''
                path = self.findShortestPath(V, E, q_init, q_goal)
                if path == None:
                    if attempt == 5:
                        break
                    n += 200
                    V, E = self.roadmapReconstruction(bounds, radius, k, n, isLazy, V, E, q_init, q_goal)
                else:
                    break

        # returns path: a list of tuples representing points along the shortest path.
        return path, V, E

    def roadmapConstruction(self, bounds, radius, k, n, isLazy):

        # Vertices in our roadmap
        V = set()
        # Edges in our roadmap
        E = set()
        '''print "start finding n free points"'''
        while len(V) < n:
            # random collision free configuration in our configuration space
            q = self.find_random_collision_free_configuration(bounds, radius)
            # V is the union of itself and the n nodes we put into our roadmap
            V.add(q)
        '''print "starting finding k nearest neighbours"'''
        for q in V:
            # k closest neighbours of configuration q, from V, sorted by distance
            Nq = self.find_k_closest_neighbours(V, q, k)
            # set of collision free edges between q and its k closest neighbours
            q_edges = self.find_collision_free_edges(q, Nq, radius, E, isLazy)
            # E is the union of itself and the newly found edges
            E = E.union(q_edges)

        # retuns V,E: the set of n Vertices and corresponding connected edges that form the probabilistic roadmap
        return V,E

    def roadmapReconstruction(self, bounds, radius, k, n, isLazy, V, E, q_init, q_goal):
        '''print "start finding n free points"'''
        V_new_points = set()
        while len(V) < n:
            # random collision free configuration in our configuration space
            q = self.find_random_collision_free_configuration(bounds, radius)
            # V is the union of itself and the n nodes we put into our roadmap
            V.add(q)
            V_new_points.add(q)
        '''print "starting finding k nearest neighbours"'''
        for q in V_new_points:
            # k closest neighbours of configuration q, from V, sorted by distance
            Nq = self.find_k_closest_neighbours(V, q, k)
            # set of collision free edges between q and its k closest neighbours
            q_edges = self.find_collision_free_edges(q, Nq, radius, E, isLazy)
            # E is the union of itself and the newly found edges
            E = E.union(q_edges)

        # k closest neighbours of start configuration q_init, from V, sorted by distance
        Nqinit = self.find_k_closest_neighbours(V, q_init, k)
        # k closest neighbours of goal configuration q_goal, from V, sorted by distance
        Nqgoal = self.find_k_closest_neighbours(V, q_goal, k)
        # V is the union of itself and Nqinit and Nqgoal
        V = V.union(Nqinit, Nqgoal)

        for qprime in Nqinit:
            #if (q_init, qprime) not in E and (qprime, q_init) not in E:
                if self.isEdgeCollisionFree(q_init, qprime, radius):
                    E.add((q_init, qprime))
                    break

        for qprime in Nqgoal:
            #if (q_goal, qprime) not in E and (qprime, q_goal) not in E:
                if self.isEdgeCollisionFree(q_goal, qprime, radius):
                    E.add((q_goal, qprime))
                    break

        return V,E

    def find_random_collision_free_configuration(self, bounds, radius):
        # continue until a point is found
        while True:
            # get a random point
            q = self.get_random_point(bounds)
            # verify if it is collision free with all obstacles and whether it is within the bounds of the problem space
            if self.isPointCollisionFree(q, radius) and not self.isOutOfBounds(q, bounds, radius) :
                return q

    def find_k_closest_neighbours(self, V, q, k):
        # convert V to a list so that we can index elements inside
        list_V = list(V)

        # initialise list to store euclidian distance between every point v in V and q
        neighbour_distances = []

        # calculate the distances and populate the list
        for v in list_V:
            neighbour_distances.append(self.euclidianDist(q, v))

        # sort the neighbours_distances list and extract their position index in V
        sorting_indexes = np.argsort(neighbour_distances)

        # extract k closest neighbours to q, by extracting the elements from V, with indexes in the 1:k+1 slots of sorting_indexes.
        # Note: We ignore the first index in sorting_indexes, because this will simply index the point q in V
        k_closest_neighbours = [list_V[sorting_indexes[i]] for i in range(1,k+1)]

        # we return a set containing the k closest neighbours, so that the union can be taken with V
        return set(k_closest_neighbours)

    # finds the set of all edges between q and its neighbours in Nq, that are collision free
    def find_collision_free_edges(self, q, Nq, radius, E, isLazy):
        # initialise set that will store all collision free edges
        edges = set()
        # iterate through all neighbours of q in Nq and find collision free edges
        for neighbour in Nq:
            if (q, neighbour) not in E and (neighbour,q) not in E:
                # add the tuple edge to edges if it is collision free. edge  = (q, neighbour)
                if isLazy:
                    edges.add((q, neighbour))
                else:
                    if self.isEdgeCollisionFree(q, neighbour, radius):
                        edges.add((q, neighbour))
        # return the set of edges so that it's union may be taken with E
        return edges

    # checks if the edge between two points in our configuration space have a collision free edge
    def isEdgeCollisionFree(self, q_init, qprime, radius):
        # generate a line between the two points using the shapely library
        line = LineString([q_init, qprime])
        # buffer the line so that it represents the physicality (shape) of the object
        expanded_line = line.buffer(radius)
        # iterate through all obstacles in the environment and verify if the buffered edge interesects any one of them
        for obstacle in self.env.obstacles:
            # terminate false as soon as intersection exists
            if expanded_line.intersects(obstacle):
                return False
        # terminate true if no intersection exists
        return True

    def findShortestPath(self, V, E, start_state, goal_state):
        graph = self.generateGraph(V,E)
        problem = GraphSearchProblem(graph, start_state, goal_state)
        return self.astar_search(problem, self.heuristic_to_goal)

    def generateGraph(self, V, E):
        graph = Graph()
        for edge in E:
            graph.add_edge(edge[0], edge[1], self.euclidianDist(edge[0], edge[1]))
        vertex_locations_dictionary = dict()
        for vertex in V:
            vertex_locations_dictionary[vertex] = vertex
        graph.set_node_positions(vertex_locations_dictionary)
        return graph

    def astar_search(self, problem, h):
        cost_func = lambda x: x.cost
        f = lambda x: cost_func(x) + h(x, problem.goal)
        return self.best_first_search(problem, f)

    def best_first_search(self, problem, f):
        queue = PriorityQueue(f=f)
        queue.append(SearchNode(problem.start))
        expanded = set([problem.start])
        max_queue = 1
        path = list()
        while queue:
            current_node = queue.pop()
            expanded.add(current_node.state)

            if(current_node.state == problem.goal):
                final_path = Path(current_node).path
                return final_path

            expanded_sn = problem.expand_node(current_node)

            for sn in expanded_sn:
                if(sn.state not in expanded):
                    if sn in queue:
                        if(sn.cost < queue[sn].cost):
                            del queue[sn]
                            queue.append(sn)
                    else:
                        queue.append(sn)
            if(len(queue) > max_queue):
                max_queue = len(queue)
        # If we get to here, no solution has been found.
        return None

    def heuristic_to_goal(self, search_node, goal_state):
        return self.euclidianDist(search_node.state, goal_state)

    # find a random point within the configuration space
    # Note: it finds an actual point and not a buffered point that represents the phsicality of our object
    def get_random_point(self, bounds):
        x = bounds[0] + random.random()*(bounds[2]-bounds[0])
        y = bounds[1] + random.random()*(bounds[3]-bounds[1])
        return (x, y)

    # check if configuration point q is collision free with all obstacles in the environment
    def isPointCollisionFree(self, q, radius):
        # buffer the point using the shapely class to check if the physical representation of q intersects any obstacles
        buffered_point = Point(q).buffer(radius)
        # verify collisions with every obstacle in the environment
        for obstacle in self.env.obstacles:
            # if q is in any obstacle, terminate false
            if obstacle.intersects(buffered_point):
                return False
        # if no obstacle contains q, then terminate true
        return True

    # check if configuration point q is out of bounds
    def isOutOfBounds(self, q, bounds, radius):
        # left boundary check
        if((q[0]-radius) < bounds[0]):
            return True
        # bottom boundary check
        if((q[1]-radius) < bounds[1]):
            return True
        # right boundary check
        if((q[0]+radius) > bounds[2]):
            return True
        # top boundary check
        if((q[1]+radius) > bounds[3]):
            return True
        # if all checks fail, then ball q wih given radius is within the boundary
        return False

    # calculate the euclidian distance between two configuration points q1, q2
    def euclidianDist(self, q1, q2):
        return math.sqrt((q2[0]-q1[0])**2 + (q2[1]-q1[1])**2)

class GraphSearchProblem(object):
    def __init__(self, graph, start, goal):
        self.graph = graph
        self.start = start
        self.goal = goal
    def test_goal(self, state):
        return self.goal == state
    def expand_node(self, search_node):
        """Return a list of SearchNodes, having the correct state, parent and updated cost."""
        current_node = search_node.state
        current_cost = search_node.cost
        outgoing_edges = self.graph.node_edges(current_node)
        expanded_sn = [SearchNode(edge.target, search_node, current_cost + edge.weight) for edge in outgoing_edges ]
        return expanded_sn


text:
import tensorflow as tf
import numpy as np
import os

batch_size = 10
img_path = "/gpfs/fs01/user/s076-844c78348e985f-04662317cedd/notebook/work/Dataset/flickr30k-images/"
try:
    files = sorted(np.array(os.listdir("/gpfs/fs01/user/s076-844c78348e985f-04662317cedd/notebook/work/Dataset/flickr30k-images/")))
    n_batch = len(files) / batch_size
except:
    pass

with open('/gpfs/global_fs01/sym_shared/YPProdSpark/user/s076-844c78348e985f-04662317cedd/notebook/work/farman-image-caption/ConvNets/inception_v4.pb', 'rb') as f:
    fileContent = f.read()

graph_def = tf.GraphDef()
graph_def.ParseFromString(fileContent)
tf.import_graph_def(graph_def)
graph = tf.get_default_graph()

input_layer = graph.get_tensor_by_name("import/InputImage:0")
output_layer = graph.get_tensor_by_name(
    "import/InceptionV4/Logits/AvgPool_1a/AvgPool:0")

'''
OLD PRE-PROCESSING MODULES : SLOW
import cv2
from PIL import Image

def old_load_image(x, new_h=299, new_w=299):
    image = Image.open(x)
    h, w = image.size
    if image.format != "PNG":
        image = np.asarray(image)/255.0
    else:
        image = np.asarray(image)/255.0
        image = image[:,:,:3]

    ##To crop or not?
    if w == h:
        resized = cv2.resize(image, (new_h,new_w))
    elif h < w:
        resized = cv2.resize(image, (int(w * float(new_h)/h), new_w))
        crop_length = int((resized.shape[1] - new_h) / 2)
        resized = resized[:,crop_length:resized.shape[1] - crop_length]
    else:
        resized = cv2.resize(image, (new_h, int(h * float(new_w) / w)))
        crop_length = int((resized.shape[0] - new_w) / 2)
        resized = resized[crop_length:resized.shape[0] - crop_length,:]

    return cv2.resize(image, (new_h, new_w))
'''

def build_prepro_graph():
    input_file = tf.placeholder(dtype=tf.string, name="InputFile")
    image_file = tf.read_file(input_file)
    jpg = tf.image.decode_jpeg(image_file, channels=3)
    png = tf.image.decode_png(image_file, channels=3)
    output_jpg = tf.image.resize_images(jpg, [299, 299]) / 255.0
    output_jpg = tf.reshape(
        output_jpg, [
            1, 299, 299, 3], name="Preprocessed_JPG")
    output_png = tf.image.resize_images(png, [299, 299]) / 255.0
    output_png = tf.reshape(
        output_png, [
            1, 299, 299, 3], name="Preprocessed_PNG")
    return input_file, output_jpg, output_png

def load_image(sess, io, image):
    if image.split('.')[-1] == "png":
        return sess.run(io[2], feed_dict={io[0]: image})
    return sess.run(io[1], feed_dict={io[0]: image})

def load_next_batch(sess, io):
    for batch_idx in range(0, len(files), batch_size):
        batch = files[batch_idx:batch_idx + batch_size]
        batch = np.array(
            map(lambda x: load_image(sess, io, img_path + x), batch))
        batch = batch.reshape((batch_size, 299, 299, 3))
        yield batch

def forward_pass(io):
    global output_layer
    with tf.Session() as sess:
        init = tf.global_variables_initializer()
        sess.run(init)
        batch_iter = load_next_batch(sess, io)
        for i in xrange(n_batch):
            batch = batch_iter.next()
            assert batch.shape == (batch_size, 299, 299, 3)
            feed_dict = {input_layer: batch}
            if i is 0:
                prob = sess.run(
                    output_layer, feed_dict=feed_dict).reshape(
                    batch_size, 1536)
            else:
                prob = np.append(
                    prob,
                    sess.run(
                        output_layer,
                        feed_dict=feed_dict).reshape(
                        batch_size,
                        1536),
                    axis=0)
            if i % 5 == 0:
                print "Progress:" + str(((i + 1) / float(n_batch) * 100)) + "%\n"
    print "Progress:" + str(((n_batch) / float(n_batch) * 100)) + "%\n"
    print
    print "Saving Features : features.npy\n"
    np.save('/gpfs/fs01/user/s076-844c78348e985f-04662317cedd/notebook/work/Dataset/features', prob)

def get_features(sess, io, img, saveencoder=False):
    global output_layer
    output_layer = tf.reshape(output_layer, [1,1536], name="Output_Features")
    image = load_image(sess, io, img)
    feed_dict = {input_layer: image}
    prob = sess.run(output_layer, feed_dict=feed_dict)

    if saveencoder:
        tensors = [n.name for n in sess.graph.as_graph_def().node]
        with open("/gpfs/fs01/user/s076-844c78348e985f-04662317cedd/notebook/work/model/Encoder/Encoder_Tensors.txt", 'w') as f:
            for t in tensors:
                f.write(t + "\n")
        saver = tf.train.Saver()
        saver.save(sess, "/gpfs/fs01/user/s076-844c78348e985f-04662317cedd/notebook/work/model/Encoder/model.ckpt")
    return prob

if __name__ == "__main__":
    print "#Images:", len(files)
    print "Extracting Features"
    io = build_prepro_graph()
    forward_pass(io)
    print "done"

ext: py

reason: import numpy


text:
import pyautogui
from PIL import ImageGrab, ImageOps
from numpy import *

def Jump():
	pyautogui.keyDown('space')

def ImageSum():
    box = [770, 380, 850, 450]
    image = ImageGrab.grab(box)
    grayImage = ImageOps.grayscale(image)
    a = array(grayImage.getcolors())
    print(a.sum())
    return a.sum()

def main():
	while True:
		if ImageSum () > 6000:
			Jump()

if __name__ == "__main__":
	main()

reason: from numpy


text:
import numpy as np
from Board import *
import random

class General():
    #
    #General was hit - mark the as not floating anymore
    #
    def hit(self, board: Board, x, y, z):
        self.floating = False
        board.General = 0
        return "General SUNK"

    #
    #Add a general randomly on the board
    #
    def add_to_board(self, board: Board):
        # Dont add more than 1 general
        if(board.General > 0):
            return

        #Get random x,y,z
        x = random.randint(1, board.x) - 1
        y = random.randint(1, board.y)-1
        z = random.randint(1, board.z)-1

        orgx = x
        orgy = y
        orgz = z

        #Loop until you can place the sub
        while True:
            if(board.board[x,y,z] == 0):
                board.board[x,y,z] = self
                board.General = 1
                self.floating = True
                return True

            #You can not add, so increment coordinate
            x = x+1
            if(x >= board.x):
                x = 0
                y = y + 1
            if(y >= board.y):
                y = 0
                z = z + 1
            if(z >= board.z):
                z = 0
            #If we tried all options - return false
            if(x == orgx and y == orgy and z == orgz ):
                return False

    #
    #__str__  Print the submarine information
    #
    def __str__(self):
        if(self.floating):
            return "**GEN**"
        else:
            return "gen"


text:
from __future__ import print_function

# Copyright (c) 2013, Roger Lew
# All rights reserved.
#
# Credits: read_daq, and _loaddata are optimized version of Chris
#          Schwarz's convert_daq.py scipy in ndaqTools

import os
import warnings

from array import array
from collections import OrderedDict, namedtuple
from fnmatch import fnmatch
from operator import attrgetter
from random import uniform
from struct import unpack

import h5py
import matplotlib.pyplot as plt
plt.rc('font', family='serif')
import numpy as np

from numpy import fromfile
from scipy import io as sio

from undaqTools.misc.base import  _size_lookup, _nptype_lookup
from undaqTools.element import Element, FrameSlice, FrameIndex, findex
from undaqTools.dynobj import DynObj
from undaqTools.misc.base import _searchsorted
from undaqTools.misc.recordtype import recordtype
from undaqTools.misc.ast import _literal_eval, _literal_repr

interpolation_wclist = \
    ['CIS_Auxiliary_Buttons',
     'SCC_Collision*',
     'SCC_Current_VDS_Frame_Count',
     'SCC_DRT_ReactionTime',
     'SCC_DynObj*',
     'SCC_Eval*',
     'SCC_False_Alarm',
     'SCC_Fog',
     'SCC_Glare_Obj*',
     'SCC_Graphics_Frame_Time',
     'SCC_HighRes_Time',
     'SCC_Init*',
     'SCC_Lane_Depart_Warn',
     'SCC_LogStreams',
     'SCC_NVES_Target_Dist',
     'SCC_StatObj*',
     'SCC_StatObj_AudioVisualState',
     'SCC_StatObj_CvedId',
     'SCC_StatObj_DataSize',
     'SCC_StatObj_Pos',
     'SCC_StatObj_SolId',
     'SCC_Time_Counter',
     'SCC_Tire_Condition',
     'SCC_TrafLight_Id',
     'SCC_TrafLight_Size',
     'SCC_TrafLight_State',
     'SCC_Trailer_Col_Det_Ob_SolId',
     'SCC_Trailer_Col_Det_Ob_Type',
     'SCC_Trailer_Col_Det_Object',
     'SCC_Trailer_Col_List_Size',
     'SCC_Trailer_Collision_Count',
     'VDS_DRV_Frame_No',
     'VDS_Frame_Count']

Info = namedtuple('Info', ['run','runinst','title','numentries',
                           'frequency','date','magic','subject',
                           'filename'])

Frame = recordtype('Frame', ['code', 'frame', 'count'])

Header = recordtype('Header', ['id','numvalues', 'name', 'units',
                               'bytes', 'rate', 'type', 'nptype',
                               'varrateflag'])

def stat(filename):
    """
    returns an Info namedtuple instance without loading data
    so the meta information can be interogatted without having
    to read the whole file

    Parameters
    ----------
    filename : string
       .daq file to investigate

    Returns
    -------
    daq : Daq instance

    Example
    -------
    >>> from pydaqTools import stat
    >>> print(stat('drive01.daq'))
    DAQ INFO
    **************************************************
    magic       7f4e3d2c
    frequency   59
    run         data reduction
    title       Nads MiniSim
    date        Mon Feb 04 12:56:17 2013
    runinst     20130204125617
    numentries  245
    subject     part12
    """
    fid = open(filename,'rb')
    read = fid.read # dots make things slow...

    # read info
    magic = read(4).encode('hex')
    title = read(120).split('\x00')[0]
    date = read(27).split('\x00')[0]
    subject = read(128).split('\x00')[0]
    run = read(128).split('\x00')[0]
    runinst = read(129).split('\x00')[0]
    numentries = unpack('i',read(4))[0]
    frequency = unpack('i',read(4))[0]
    filename = filename

    fid.close()

    return \
        Info(run = run,
             runinst = runinst,
             title = title,
             numentries = numentries,
             frequency = frequency,
             date = date,
             magic = magic,
             subject = subject,
             filename = filename)

class Daq(dict):
    def __init__(self):
        """Abstraction of NADS .daq data"""

        # namedtuple containing drive relevant metadata
        self.info = Info(run = '',
                         runinst = '',
                         title = '',
                         numentries = 0,
                         frequency = 0,
                         date = '',
                         magic = '',
                         subject = '',
                         filename = '')

        # recordtype of frame information
        self.frame = Frame(code = array('i'),
                           frame = array('i'),
                           count = array('i'))

        # recordtype of element header infomation
        self._header = None

        # other attributes
        self.data = {}               # dictionary to store elements
        self.elemlist = None         # specifies elements to load
        self.f0 = None               # first Frame in Daq. Probably not 0
        self.fend = None             # last Frame in Daq
        self.cursor = 0              # byte where the data frames begin in daq
        self.dynobjs = OrderedDict() # object to hold dynamic objects
        self.etc = {}

        dict.__init__(self)

    def load_elemlist_fromfile(self, filename):
        with open(filename,'r') as f:
            self.elemlist = [s.rstrip() for s in f.readlines()]

    def read_daq(self, filename, elemlist=None,
                 loaddata=True, process_dynobjs=True,
                 interpolate_missing_frames=True):
        """
        read_daq(filename[, elemlist=None]
                 [, loaddata=True][, process_dynobjs=True])

        Reads a .daq file into object

        read <==> read_daq

        Parameters
        ----------
        filename : string
            path to .daq file

        elemlist : None or list_like
             None -> load all elements
             list_like -> load names that match list

        loaddata : bool
             True -> load data
             False -> don't load data

        process_dynobjs : bool
             True -> process dynobjs and put them in self.dynobjs
             False -> don't load dynobjs
        """

        _header = \
            Header(id = None,
                   numvalues = array('i'),
                   name = [],
                   units = [],
                   rate = array('i'),
                   type = array('c'),
                   nptype = [],
                   varrateflag = array('i'),
                   bytes = array('i'))

        if elemlist is not None:
            self.elemlist = elemlist

        fid = open(filename,'rb')
        read = fid.read # dots make things slow...

        # read info
        magic = read(4).encode('hex')
        title = read(120).split('\x00')[0]
        date = read(27).split('\x00')[0]
        subject = read(128).split('\x00')[0]
        run = read(128).split('\x00')[0]
        runinst = read(129).split('\x00')[0]
        numentries = unpack('i',read(4))[0]
        frequency = unpack('i',read(4))[0]
        filename = filename

        self.info = \
            Info(run = run,
                 runinst = runinst,
                 title = title,
                 numentries = numentries,
                 frequency = frequency,
                 date = date,
                 magic = magic,
                 subject = subject,
                 filename = filename)

        # read header
        _header.id = range(self.info.numentries)

        for i in _header.id:
            # this use to be append_header
            # function call overhead is kind of high with
            # python so it is moved here for performance

            _header.numvalues.append(unpack('i',read(4))[0])
            _header.name.append(read(36).split('\x00')[0])
            _header.units.append(read(16).split('\x00')[0])

            _header.rate.append(unpack('h', read(2))[0])
            read(2) # seek is slow...
            _header.type.append(chr(unpack('i', read(4))[0]))
            if _header.type[-1] == 's':
                _header.type[-1] = 'h'
            _header.varrateflag.append(unpack('B', read(1))[0])
            read(3) # seek is slow...

        _header.bytes = [_size_lookup[t] for t in _header.type]
        _header.nptype = [_nptype_lookup[t] for t in _header.type]

        # remember where we are so we know where to start
        # reading in the data
        self.cursor = fid.tell()
        self._header = _header

        fid.close()

        if loaddata:
            self._loaddata()
            self._unwrap_lane_deviation()

        if loaddata and process_dynobjs:
            self._process_dynobjs()

        # interpolate if frames are missing
        # we want to do this after loading dyn_objects so that the
        # SCC_DynObjs* stay aligned with the rest of the data.
        # DynObj.process detects whether frames are missing and
        # interpolates the dynamic object Elements so in the end
        # everything gets realigned.

        if len(self) == 0:
            return # no data loaded

        missing_frames = (self.fend - self.f0 + 1) - len(self.frame.frame)
        if missing_frames != 0:
            msg = "Missing %i frames."%missing_frames
            if interpolate_missing_frames:
                msg += " (interpolated missing frames)"
                self._interpolate_missing_frames()
            warnings.warn(msg, RuntimeWarning)

    # create alias read for read_daq
    read = read_daq

    def _loaddata(self):
        """
        loads data from .daq file
        """
        _header = self._header
        elemlist = self.elemlist
        frame = self.frame

        def nans(numitems):
            return [np.nan for i in xrange(numitems)]

        _elemid_lookup = dict(zip(_header.name, _header.id))

        # Data gets unpacked to a temporary dict tmpdata before building
        # Element objects
        #
        # We will unpack 1d elements to an array object and then
        # change them to 2d ndarrays after the data is loaded.
        #
        # Elements with more than 1 dimension will get unpacked
        # as lists of arrays and transformed to ndarrays after
        # the data is loaded
        #
        # The array.array objects are typed unlike regular lists.
        # This improves memory usage. The array.array objects also
        # have highly efficient appending unlike ndarrays

        # initialize temporary dict
        tmpdata = {}
        for name, typ, varrate, numvalues in \
            zip(_header.name, _header.type,
                _header.varrateflag, _header.numvalues):
            if elemlist is None or any(fnmatch(name, wc) for wc in elemlist):
                if not varrate and numvalues == 1:
                    tmpdata[name] = array(typ)
                else:
                    tmpdata[name] = []

        # Intialize Frame arrays for CSSDC measures
        for name in tmpdata.keys():
            if _header.rate[_elemid_lookup[name]] != 1:
                tmpdata[name+'_Frames'] = array('i')

        fid = open(self.info.filename,'rb')
        read = fid.read # dots make things slow...
        read(self.cursor) # move to where the data begins

        # The daq files have a frame for every sample. Each frame
        # contains the variables states at the time the frame was
        # saved.
        #
        # we want to create arrays that can be indexed to decide
        # whether the variable needs stored. This way we don't
        # have to look through the elemlist for every variable on
        # every frame
        if elemlist is None: # elemlist empty
            mask = [True for i in xrange(len(_header.name))]
        else:
            mask = []
            for name in _header.name:
                mask.append(any(fnmatch(name, wc) for wc in elemlist))

        bombed = False
        # Python quirk no. 372837462: while 1 is faster than while True
        while 1:
            try:
                frame.code.append(unpack('i',read(4))[0])
                if frame.code[-1] == -2:
                    break

                frame.frame.append(unpack('i', read(4))[0])
                frame.count.append(unpack('i', read(4))[0])

                # xrange (2.7, range in 3) is faster than range
                for j in xrange(frame.count[-1]):
                    i = unpack('i', read(4))[0]

                    if _header.varrateflag[i]:
                        numitems = unpack('i', read(4))[0]
                    else:
                        numitems = _header.numvalues[i]

                    size = _header.bytes[i]

                    if mask[i]:
                        name = _header.name[i]

                        if numitems == 1:
                            typ = _header.type[i]
                            tmpdata[name].append(unpack(typ,read(size))[0])

                        else: # numitems > 1
                            typ = _header.nptype[i]
                            tmpdata[name].append(fromfile(fid,typ,numitems))

                        if _header.rate[i] != 1:
                            tmpdata[name+'_Frames'].append(frame.frame[-1])

                    else: # we don't need to read this element
                        read(numitems*size) # seek is slow...
                                            # reduce calls to read
            except:
                msg = 'Failed loading file on frame %i.'%frame.frame[-1]
                msg += ' (stopped reading file)'
                warnings.warn(msg, RuntimeWarning)
                bombed = True
                break

        fid.close()

        # We made it through the daq file.
        # Now it is time to do some bookkeeping.
        frame.frame = np.array(frame.frame)
        self.f0 = f0 = frame.frame[0]
        self.fend = fend = frame.frame[-1]

        # If we bombed unpacking a frame we need to make sure
        # everything is aligned before we intialize Elements
        if bombed:
            # we will strip off the last frame just to make sure
            # everything is kosher
            n = min(len(frame.code),len(frame.frame),len(frame.count))-1
            frame.code = frame.code[:n]
            frame.frame = frame.frame[:n]
            frame.count = frame.count[:n]

            for name in tmpdata:
                if name.endswith('_Frames'):
                    continue

                val = tmpdata[name]
                i = _elemid_lookup[name]
                if _header.rate[i] != 1:
                    # CSSDC measures should always be okay. If they have 1
                    # value they use read, if they have multiple values they
                    # use fromfile. In either case if the unpacking fails
                    # nothing will get appended to tmpdata and it will exit
                    # before appending to the cooresponding _Frames array.
                    #
                    # At least that is what my mental interpreter thinks.
                    # we will assert just to be on the safe side.
                    assert len(tmpdata[name]) == len(tmpdata[name+'_Frames'])

                else:
                    tmpdata[name] = np.array(tmpdata[name][:n], ndmin=2)

                    # gets transposed in next for loop if numvalues > 1
                    if _header.numvalues[i] > 1:
                        assert tmpdata[name].shape[0] == frame.frame.shape[0]
                    else:
                        assert tmpdata[name].shape[1] == frame.frame.shape[0]

        # cast as Element objects
        # 'varrateflag' variables remain lists of lists
        #
        # There are obvious more compact ways to write this but I'm
        # paranoid about reference counting and garbage collection not
        # functioning properly
        for name, i, rate in zip(_header.name, _header.id, _header.rate):
            if elemlist is not None:
                if not any(fnmatch(name, wc) for wc in elemlist):
                    continue

            # transpose Element with more than 1 row
            if _header.numvalues[i] > 1:
                tmpdata[name] = np.array(tmpdata[name]).transpose()

            if rate != 1:
                self[name] = \
                    Element(tmpdata[name],
                            tmpdata[name+'_Frames'],
                            rate=_header.rate[i],
                            name=_header.name[i],
                            dtype=_header.type[i],
                            varrateflag=_header.varrateflag[i],
                            elemid=_header.id[i],
                            units=_header.units[i])

                del tmpdata[name+'_Frames']
            else:
                self[name] = \
                    Element(tmpdata[name],
                            frame.frame[:],
                            rate=_header.rate[i],
                            name=_header.name[i],
                            dtype=_header.type[i],
                            varrateflag=_header.varrateflag[i],
                            elemid=_header.id[i],
                            units=_header.units[i])

                # delete tmpdata arrays as we go to save memory
                del tmpdata[name]

        del _header, self._header
        self._header = None

    def _interpolate_missing_frames(self):
        """
        interpolates over missing frames for non-CSSDC measures

        Excludes elements that wildcard match:

            CIS_Auxiliary_Buttons
            SCC_Collision*
            SCC_Current_VDS_Frame_Count
            SCC_DRT_ReactionTime
            SCC_DynObj*
            SCC_Eval*
            SCC_False_Alarm
            SCC_Fog
            SCC_Glare_Obj*
            SCC_Graphics_Frame_Time
            SCC_HighRes_Time
            SCC_Init*
            SCC_Lane_Depart_Warn
            SCC_LogStreams
            SCC_NVES_Target_Dist
            SCC_StatObj*
            SCC_StatObj_AudioVisualState
            SCC_StatObj_CvedId
            SCC_StatObj_DataSize
            SCC_StatObj_Pos
            SCC_StatObj_SolId
            SCC_Time_Counter
            SCC_Tire_Condition
            SCC_TrafLight_Id
            SCC_TrafLight_Size
            SCC_TrafLight_State
            SCC_Trailer_Col_Det_Ob_SolId
            SCC_Trailer_Col_Det_Ob_Type
            SCC_Trailer_Col_Det_Object
            SCC_Trailer_Col_List_Size
            SCC_Trailer_Collision_Count
            VDS_DRV_Frame_No
            VDS_Frame_Count
        """
        global interpolation_wclist
        f0, fend = self.f0, self.fend
        old_frames = self.frame.frame[:]
        new_frames = np.linspace(f0, fend, fend - f0 + 1)

        # np.interp can only handle 1-d arrays
        for elem in self.values():
            if elem.isCSSDC():
                continue

            if any(fnmatch(elem.name, wc) for wc in interpolation_wclist):
                # The 'SCC_DynObj*' cells are non-CSSDC but contain
                # categorical data so we don't want to interpolate it.
                #
                # All we really have to do is tell Daq to treat it as
                # a non-CSSDC measures. When it gets exported the
                # frames will go with it as SCC_DynObj*_Frames
                self[elem.name].rate = -2
                continue

            x = [np.interp(new_frames, old_frames, elem[j,:])
                 for j in xrange(elem.numvalues)]

            self[elem.name] = \
                Element(x, new_frames,
                        rate=elem.rate,
                        varrateflag=elem.varrateflag,
                        elemid=elem.id,
                        units=elem.units)

        self.frame.frame = new_frames

        # did we fix the problem?
        assert fend - f0 + 1 == len(new_frames)

    def _process_dynobjs(self):
        """
        figure out what rows and frames relate to particular
        dynamic objects so that they can be unpacked by
        DynObj.process
        """

        if 'SCC_DynObj_CvedId'not in self:
            msg = "Need 'SCC_DynObj*' to build dynobjs (did not build dynobjs)"
            warnings.warn(msg, RuntimeWarning)
            return None

        frame = self.frame
        cvedid_array = self['SCC_DynObj_CvedId']

        # identify the cars in the datafile and
        # identify what frames they are defined
        frame_indiceses = OrderedDict()
        row_indiceses = OrderedDict()

        # loop through frames, skipping first frame
        for i in xrange(1, len(frame.frame)):

            cveds_in_frame = cvedid_array[:,i].flatten()

            for j, cved in enumerate(cveds_in_frame):
                if cved != 0:
                    try:
                        # if cved not in dict it will raise a
                        # key error
                        frame_indiceses[cved].append(i)
                        row_indiceses[cved].append(j)
                    except:
                        frame_indiceses[cved] = array('i', [i])
                        row_indiceses[cved] = array('i', [j])

        self.dynobjs = OrderedDict()
        for cvedId in frame_indiceses.keys():
            do = DynObj()
            do.process(cvedId,
                       np.array(frame_indiceses[cvedId]),
                       np.array(row_indiceses[cvedId]),
                       self)
            self.dynobjs[do.name] = do

        del frame_indiceses
        del row_indiceses

    def __setitem__(self, name, elem):
        if not isinstance(elem, Element):
            raise(TypeError, 'Value must be Element')

        if name is None and elem.name is not None:
            name = elem.name

        if elem.name is None and name is not None:
            elem.name = name

        if name != elem.name:
            msg = "key:'%s' does not equal elem.name'%s'\nCannot set value."\
                  %(name, elem.name)
            raise(KeyError, msg)

        if name is None:
            msg = "name cannot be None"
            raise(KeyError, msg)

        dict.__setitem__(self, name, elem)

    def _unwrap_lane_deviation(self, lanewidth=12.):
        """
        Uses np.unwrap to create a continuous lane deviation measure

        measure is registered in self as 'SCC_Spline_Lane_Deviation_Fixed'

        Has not been extensively tested or evaluated.

        Parameters
        ----------
        lanewidth : float
            assumed lanewidth in feet
        """
        if 'SCC_Spline_Lane_Deviation'not in self:
            msg = "Need 'SCC_Spline_Lane_Deviation' to fix lane deviation"
            msg += ' (did not fix lane deviation)'
            warnings.warn(msg, RuntimeWarning)
            return None

        halfwidth = lanewidth/2.
        lane_dev = self['SCC_Spline_Lane_Deviation'][1,:]

        fixed_lane_dev = lane_dev*(np.pi/halfwidth) # convert to radians
        fixed_lane_dev = np.unwrap(fixed_lane_dev)  # unwrap
        fixed_lane_dev *= (halfwidth/np.pi)         # convert back to feet

        self['SCC_Spline_Lane_Deviation_Fixed'] = \
            Element(fixed_lane_dev,
                    self.frame.frame[:],
                    name=lane_dev.name+'_Fixed',
                    elemid=lane_dev.id,
                    rate=lane_dev.rate,
                    varrateflag=lane_dev.varrateflag,
                    units=lane_dev.units,
                    dtype=lane_dev.type)

    def _rebuild_header(self):
        """
        To write the .mat files and .hd5 files we need the header arrays.
        But inside the Daq object it is much simplier to access the element
        meta data through the Element objects. The HDF5 files could have
        a format that more closely approximates how the Daq is organized but
        it would not be as efficient storage-wise because every element would
        require have a frames array.

        If users manually add an Element to the data dict the header would
        also need to track those changes and that introduces a host of possible
        errors. It seems much simplier to just rebuild the header arrays when
        they are needed.
        """
        global _nptype_lookup, _size_lookup
        _header = \
            Header(id = array('i'),
                   numvalues = array('i'),
                   name = [],
                   units = [],
                   rate = array('i'),
                   type = array('c'),
                   nptype = [],
                   varrateflag = array('i'),
                   bytes = array('i'))

        i = 0
        for elem in sorted(self.values(), key=attrgetter('name')):
            _header.id.append(i)
            _header.numvalues.append(elem.numvalues)
            _header.name.append(elem.name)
            _header.units.append(elem.units)
            _header.rate.append(elem.rate)
            _header.type.append(elem.type)
            _header.nptype.append(elem.nptype)
            _header.varrateflag.append(elem.varrateflag)
            _header.bytes.append(elem.bytes)
            i += 1

        return _header

    def read_hd5(self, filename, f0=None, fend=None):
        """
        read_hd5(filename[, f0=None][, fend=None])

        writes Daq object to HDF5 container

        f0 and fend specify frame range to read.
        Will export daq.dynobjs to HDF5 container

        Parameters
        ----------
        filename : string
            file to read

        f0 : None or int
            None -> read from beginning of file
            int -> read from this frame

        fend : None or int
            None -> read to end of file
            int -> read to this frame

        Return
        ------
        None
        """
        if filename.endswith('daq'):
            # Normal Python mantra usually goes against this sort of
            # of checking but when you make this error the traceback
            # is somewhat obtuse. It can be difficult to tell whether
            # it is in the wrong format or something else is wrong.
            raise ValueError('You are trying to open a .daq as .hd5')

        root = h5py.File(filename, 'r')

        # info
        self.info = \
            Info(run = root['info'].attrs['run'],
                 runinst = root['info'].attrs['runinst'],
                 title = root['info'].attrs['title'],
                 numentries = root['info'].attrs['numentries'],
                 frequency = root['info'].attrs['frequency'],
                 date = root['info'].attrs['date'],
                 magic = root['info'].attrs['magic'],
                 subject = root['info'].attrs['subject'],
                 filename = root['info'].attrs['filename'])

        # header
        # The [:] unpacks the data from a h5py.dataset.Dataset
        # object to a numpy.ndarray object
        try:
            _header = \
                Header(id = root['header/id'][:],
                       numvalues = root['header/numvalues'][:],
                       name = root['header/name'][:],
                       units = root['header/units'][:],
                       rate = root['header/rate'][:],
                       type = root['header/type'][:],
                       nptype = root['header/nptype'][:],
                       varrateflag = root['header/varrateflag'][:],
                       bytes = root['header/bytes'][:])
        except:
            _header = \
                Header(id = array('i'),
                       numvalues = array('i'),
                       name = [],
                       units = [],
                       rate = array('i'),
                       type = array('c'),
                       nptype = [],
                       varrateflag = array('i'),
                       bytes = array('i'))

        # Find the indices cooresponding to the first and last
        # frames requested. We can use these indices to
        # slice out all the non-CSSDC elements.
        #
        # For the CSSDC elements we will have to find the
        # appropriate indices as we go.
        i0, iend = None, None

        try:
            all_frames = root['frame/frame'][:]
        except:
            all_frames = None

        if all_frames is not None:
            if f0 is not None or fend is not None:
                if f0 is not None:
                    i0 = _searchsorted(all_frames, f0)
                if fend is not None:
                    iend = _searchsorted(all_frames, fend)
                    if iend < len(all_frames):
                        iend += 1

        indx = slice(i0,iend)

        # frame
        try:
            self.frame = \
                Frame(code = root['frame/code'][indx],
                      frame = root['frame/frame'][indx],
                      count = root['frame/count'][indx])
        except:
            self.frame = \
                Frame(code = array('i'),
                      frame = array('i'),
                      count = array('i'))

        # elemlist
        try:
            # Fails if slice is zero-length
            self.elemlist = root['elemlist'][:]
        except:
            self.elemlist = None

        # data
        # Procedure is similar to read_daq. Data is unpacked
        # to tmpdata dict and then the Elements are instantiated.
        _elemid_lookup = dict(zip(_header.name, _header.id))

        tmpdata = {}
        for k, v in root['data'].iteritems():

            if self.elemlist is not None:
               if not any(fnmatch(k, wc) for wc in self.elemlist):
                   continue

            i = _elemid_lookup[k.replace('_Frames','')]

            if _header.rate[i] == 1:
                tmpdata[k] = v[:,i0:iend]

            else: #CSSDC measure
                if len(v.shape) == 1:
                    v = np.array(v, ndmin=2) # _Frames

                # Need to find indices
                _i0 = 0
                _iend= v.shape[1]

                if f0 is not None or fend is not None:
                    _name = k.replace('_Frames','')
                    _all_frames = root['data/%s_Frames'%_name][:].flatten()

                    if f0 is not None:
                        _i0 = _searchsorted(_all_frames, f0)
                    if fend is not None:
                        _iend = _searchsorted(_all_frames, fend)
                        if _iend < len(_all_frames):
                            _iend += 1

                # Now we can slice the data
                tmpdata[k] = v[:,_i0:_iend]

        # hdf5 doesn't have a None type (or atleast, I don't know how
        # to use it) so None is stored as an empty string in the hdf5 file
        self.f0 = (root.attrs['f0'], None)[root.attrs['f0'] == '']
        self.fend = (root.attrs['fend'], None)[root.attrs['fend'] == '']
        self.cursor = root.attrs['cursor']

        # read dynobjs
        self.dynobjs = OrderedDict()
        for (name, dynobj) in root['dynobjs'].iteritems():
            do = DynObj()
            do.read_hd5(root=dynobj)
            self.dynobjs[name] = do

        # read etc dict
        self.etc = {}
        for (name, obj) in root['etc'].attrs.iteritems():
            self.etc[name] = _literal_eval(obj)

        root.close()

        # cast as Element objects
        # 'varrateflag' variables remain lists of lists
        #
        # There are obvious more compact ways to write this but I'm
        # paranoid about reference counting and garbage collection not
        # functioning properly
        for name, i, rate in zip(_header.name, _header.id, _header.rate):
            if rate != 1:
                self[name] = \
                    Element(tmpdata[name],
                            tmpdata[name+'_Frames'].flatten(),
                            rate=_header.rate[i],
                            name=_header.name[i],
                            dtype=_header.type[i],
                            varrateflag=_header.varrateflag[i],
                            elemid=_header.id[i],
                            units=_header.units[i])

                del tmpdata[name+'_Frames']
            else:
                self[name] = \
                    Element(tmpdata[name],
                            self.frame.frame[:],
                            rate=_header.rate[i],
                            name=_header.name[i],
                            dtype=_header.type[i],
                            varrateflag=_header.varrateflag[i],
                            elemid=_header.id[i],
                            units=_header.units[i])

                # delete tmpdata arrays as we go to save memory
                del tmpdata[name]

        del _header

    def write_hd5(self, filename=None):
        """
        write_hd5(filename=None)

        writes Daq object to HDF5 container

        Parameters
        ----------
        filename : None or string (optional)
            None : file written to daq.filename.replace('.daq', '.hdf5')
            string : specify output file (will overwrite)

        Return
        ------
        None
        """
        info = self.info
        frame = self.frame
        _header = self._rebuild_header()
        str_type = h5py.new_vlen(str)

        if filename is None:
            filename = info.filename[:-4] + '.hdf5'

        if filename.endswith('.daq'):
            # yes, I did this accidently
            msg = 'Writing Daq in HDF5 with a .daq extension is not allowed'
            raise ValueError(msg)

        root = h5py.File(filename, 'w')

        # info
        root.create_group('info')
        root['info'].attrs['run'] = info.run
        root['info'].attrs['runinst'] = info.runinst
        root['info'].attrs['title'] = info.title
        root['info'].attrs['numentries'] = info.numentries
        root['info'].attrs['frequency'] = info.frequency
        root['info'].attrs['date'] = info.date
        root['info'].attrs['magic'] = info.magic
        root['info'].attrs['subject'] = info.subject
        root['info'].attrs['filename'] = info.filename

        # frame
        root.create_group('frame')
        root['frame'].create_dataset('frame', data=frame.frame)
        root['frame'].create_dataset('count', data=frame.count)
        root['frame'].create_dataset('code', data=frame.code)

        # header
        root.create_group('header')
        root['header'].create_dataset('id', data=_header.id)
        root['header'].create_dataset('numvalues', data=_header.numvalues)
        root['header'].create_dataset('name', data=_header.name)
        root['header'].create_dataset('units', data=_header.units)
        root['header'].create_dataset('rate', data=_header.rate)
        root['header'].create_dataset('type', data=_header.type)
        root['header'].create_dataset('varrateflag', data=_header.varrateflag)
        root['header'].create_dataset('bytes', data=_header.bytes)

        nptype_asstr = np.array(_header.nptype, dtype=str_type)
        root['header'].create_dataset('nptype', data=nptype_asstr)

        # data
        root.create_group('data')
        for name, elem in self.items():
            root['data'].create_dataset(name, data=elem.toarray())

            if elem.isCSSDC():
                root['data'].create_dataset(name+'_Frames',
                                            data=elem.frames)
        # misc
        root.create_dataset('elemlist',
                            data=np.array(self.elemlist, dtype=str_type))
        root.attrs['f0'] = (self.f0, '')[self.f0 is None]
        root.attrs['fend'] = (self.fend, '')[self.fend is None]
        root.attrs['cursor'] = self.cursor

        # dynobjects
        root.create_group('dynobjs')
        for (name, do) in self.dynobjs.items():
            root.create_group('dynobjs/%s'%name)
            do.write_hd5(root=root['dynobjs/%s'%name])

        # etc
        root.create_group('etc')
        for (name, obj) in self.etc.items():
            s = _literal_repr(obj)

            try:
                _literal_eval(s)
            except:
                root.close()
                msg = "'%s' is not a Python literal, "\
                      "FrameSlice, or FrameIndex"%str(obj)
                raise ValueError(msg)
            root['etc'].attrs[name] = s

        root.close()

    def write_mat(self, filename=None, outpath=None):
        """
        write_mat([filename=None][, outpath=None])

        writes Daq object to .mat file using scipy.io.savmat

        daq.dynobjs is not exported

        This is intended to make convert_daq backwards compatible.
        pydaqTools users should use write_hdf5. The .mat files cannot
        be used to instantiate Daq objects.

        Parameters
        ----------
        filename : None or string (optional)
            None : file written to daq.filename.replace('.daq', '.hdf5')
            string : specify output file (will overwrite)

        outpath : None or string (optional)
            string : directory to path

        Return
        ------
        None
        """

        if filename is None:
            filename = os.path.splitext(self.info.filename)[0]+'.mat'

        if outpath is not None:
            filename = os.path.join(outpath,filename)

        # data
        data = {}
        for name in self:
            data[name] = np.array(self[name], order='F').transpose()

            if self[name].isCSSDC():
                data[name+'_Frames'] = np.array(self[name].frames, order='F')
                data[name+'_Frames'] = data[name+'_Frames'].transpose()

        data['Frames'] = np.array(self.frame.frame, ndmin=2, order='F')
        data['Frames'] = data['Frames'].transpose()

        # header
        header = self._rebuild_header()._asdict()
        del header['nptype']
        header['type'] = [(t,'s')[t=='h'] for t in header['type']]

        sio.savemat(filename,
                    {'daqInfo':dict(self.info._asdict()),
                    'elemInfo':header,
                    'elemData':data,
                  'elemFrames':dict(self.frame._asdict())},
                    long_field_names=True, oned_as='column')

        del data
        del header

    def plot_ts(self, elem_pars, xindx=None):
        """
        plot_ts(self, elem_pars[, xindx=None])

        Returns a time series figure containing elements specified by
        elem_pars

        CSSDC measures yield step plots.

        Parameters
        ----------
        elem_pars : list of tuples
            each tuple should contain
                name : str
                    element to plot
                yindx : None, int, or slice object
                    specifies the values of element to plot
                    None -> slice(None, None, None)

        xindx : None, slice object, FrameSlice object
            specifies the frame range to plot
            None -> slice(None, None, None)

        Returns
        -------
        fig : matplotlib.figure.Figure
            the plot
        """
        # initialize variables and figure

        fs = float(self.info.frequency)
        num_subplots = len(elem_pars)
        fig, axs = plt.figure(figsize=(16, 2*len(elem_pars))), []
        fig.subplots_adjust(left=.07, right=.98, bottom=.03,
                            top=.97, hspace=.25)

        if xindx is None:
            xindx = slice(None)

        # loop through requested elements and plot
        for i, (name, yindx) in enumerate(elem_pars):

            # make subplot
            axs.append(fig.add_subplot(num_subplots, 1, i+1))

            if yindx is None:
                yindx = slice(None)

            # figure out data
            elem = self[name]
            y = elem[yindx, xindx]
            x = y.frames.flatten()/fs

            if len(y.shape) == 1:
                if elem.isCSSDC():
                    axs[-1].step(x, y, where='pre')
                else:
                    axs[-1].plot(x, y)
            else:
                for j in xrange(y.shape[0]):
                    if elem.isCSSDC():
                        axs[-1].step(x, y[j,:], where='pre')
                    else:
                        axs[-1].plot(x, y[j,:])

            # set xlim
            axs[-1].set_xlim([x[0], x[-1]])

            # put a title on the subplot
            def _frmttr(slice_obj):
                if slice_obj == slice(None):
                    return ':'

                return str(slice_obj).replace('FrameSlice(', '')\
                                     .replace('slice(', '')\
                                     .replace('start=None', '')\
                                     .replace('stop=None', '')\
                                     .replace('step=', '')\
                                     .replace('start=', 'f')\
                                     .replace('stop=', 'f')\
                                     .replace(', None)', '')\
                                     .replace(', ', ':')\
                                     .replace(')', '')

            title = '%s[%s, %s]'%(name, _frmttr(yindx), _frmttr(xindx))
            axs[-1].set_title(title)

        # remove xticklabels from all but the last subplot
        for i in xrange(num_subplots-1):
            axs[i].set_xticklabels([])

        return fig

    def plot_dynobjs(self, dynobj_wc='*'):
        """
        plot_dynobjs(self[, dynobj_wc='*'])

        make top-down plot of dynamic object paths

        Parameters
        ----------
        dynobj_wc : string
            wildcard pattern to apply to dynamic object names using fnmatch
            (case normalized, shell-style matching)

        Returns
        -------
        fig : matplotlib.figure.Figure
            the plot
        """

        # make top down plot
        fig = plt.figure(figsize=(12, 12))
        fig.subplots_adjust(left=.02, right=.98, bottom=.02, top=.98)
        ax = fig.add_subplot(111)
        ax.axis('equal')

        colors = 'bgrcm'
        line_styles = ['-','-.',':']
        ados = [do for do in self.dynobjs.values() \
                if fnmatch(do.name, dynobj_wc)]
        for i, do in enumerate(ados):
            ax.plot(do.pos[0,:] + uniform(-600,600),
                    do.pos[1,:] + uniform(-600,600),
                    label=do.name, alpha=1.,
                    c=colors[i%5],
                    ls=line_styles[i%3])

        # plot the external driver
        ax.plot(self['VDS_Chassis_CG_Position'][0,:],
                self['VDS_Chassis_CG_Position'][1,:],
                c='k', label=self.info.subject)

##        xlim=pylab.xlim()
##        pylab.xlim([xlim[0], xlim[1]*1.9])
##        pylab.legend(ncol=2)

        return fig

    def match_keys(self, wc):
        """
        Returns a list of keys (element names) that match
        the Unix-style wildcard.

        Parameters
        ----------
        wc : string
            wildcard apply to each match

        Returns
        -------
        matches : list
            names (keys to daq)

        Notes
        -----
        keys and wc both case normalized

        Patterns are Unix shell style::

            *       matches everything
            ?       matches any single character
            [seq]   matches any character in seq
            [!seq]  matches any char not in seq

        """
        return [name for name in self if fnmatch(name, wc)]

    def keys_summary(self, wclist=None):
        """
        Print a table summary of the Elements in daq

        Parameters
        ----------
        wclist : list of strings
            list of wildcards to match to each key in daq
        """

        if wclist is None:
            wclist = ['*']

        if len(self) == 0:
            print('daq is empty')
            return

        names =  []
        for name in self:
            if any(fnmatch(name, wc) for wc in wclist):
                names.append(name)

        ncol = max(len(n) for n in names) + 1
        missing = self.fend - self.f0 + 1 - len(self.frame.frame)

        for name in sorted(names):
            elem = self[name]
            if elem.isCSSDC():
                continue

            print(name.ljust(ncol),
                  ('     ','X    ')[elem.isCSSDC()],
                  str(elem.dtype).ljust(10),
                  elem[0,:10].toarray())

max_stars_repo_name: rogerlew/undaqTools

max_issues_repo_name: rogerlew/undaqTools

max_forks_repo_name: rogerlew/undaqTools

reason: import numpy,from numpy,from scipy


text:
"""Script to disaggregate national data into regional data
"""
import os
import numpy as np
from energy_demand.assumptions import assumptions
from energy_demand.read_write import data_loader
from energy_demand.profiles import hdd_cdd

'''
============================================
MEthod to derive GVA/POP SERVICE FLOOR AREAS
============================================

1. Step
Get correlation between regional GVA and (regional floor area/reg pop) of every sector of base year
-- Get this correlation for every region and build national correlation

2. Step
Calculate future regional floor area demand based on GVA and pop projection

'''

def disaggregate_base_demand(data):
    """This function disaggregates fuel demand based on region specific parameters
    for the base year

    The residential, service and industry demand is disaggregated according to
    different factors

    Parameters
    ----------
    data : dict
        Contains all data not provided externally

    Returns
    -------
    data : dict
    """
    def sum_fuels_before(fuel):
        """Inner function for testing purposes - sum fuel"""
        tot = 0
        for i in fuel:
            tot += np.sum(fuel[i])
        return tot

    def sum_fuels_after(reg_fuel):
        """Inner function for testing purposes - sum fuel"""
        tot = 0
        for reg in reg_fuel:
            for enduse in reg_fuel[reg]:
                tot += np.sum(reg_fuel[reg][enduse])
        return tot

    # Disaggregate residential submodel data
    data['rs_fueldata_disagg'] = rs_disaggregate(data, data['rs_fuel_raw_data_enduses'])

    # Disaggregate service submodel data
    data['ss_fueldata_disagg'] = ss_disaggregate(data, data['ss_fuel_raw_data_enduses'])

    # Disaggregate industry submodel data
    data['is_fueldata_disagg'] = is_disaggregate(data, data['is_fuel_raw_data_enduses'])

    # Disaggregate transportation sector
    data['ts_fueldata_disagg'] = scrap_ts_disaggregate(data, data['ts_fuel_raw_data_enduses'])

    # Check if total fuel is the same before and after aggregation
    test_sum_before = sum_fuels_before(data['rs_fuel_raw_data_enduses'])
    test_sum_after = sum_fuels_after(data['rs_fueldata_disagg'])
    np.testing.assert_almost_equal(test_sum_before, test_sum_after, decimal=2, err_msg="")

    return data

def ss_disaggregate(data, raw_fuel_sectors_enduses):
    """Disaggregate fuel for service submodel (per enduse and sector)
    """
    print("... disaggregate service demand")
    ss_fueldata_disagg = {}

    # ---------------------------------------
    # Calculate heating degree days for regions
    # ---------------------------------------
    ss_hdd_individ_region = hdd_cdd.get_hdd_country(data['lu_reg'], data, 'ss_t_base_heating')

    ss_cdd_individ_region = hdd_cdd.get_cdd_country(data['lu_reg'], data, 'ss_t_base_cooling')

    # ---------------------------------------
    # Overall disaggregation factors per enduse and sector
    # ---------------------------------------
    # Total floor area for every enduse per sector
    national_floorarea_by_sector = {}
    for sector in data['ss_sectors']:
        national_floorarea_by_sector[sector] = 0
        for region in data['lu_reg']:
            national_floorarea_by_sector[sector] += data['ss_sector_floor_area_by'][region][sector]

    f_ss_catering = {}
    f_ss_computing = {}
    f_ss_cooling_ventilation = {}
    f_ss_water_heating = {}
    f_ss_space_heating = {}
    f_ss_lighting = {}
    f_ss_other_electricity = {}
    f_ss_other_gas = {}

    for sector in data['all_sectors']:
        f_ss_catering[sector] = 0
        f_ss_computing[sector] = 0
        f_ss_cooling_ventilation[sector] = 0
        f_ss_water_heating[sector] = 0
        f_ss_space_heating[sector] = 0
        f_ss_lighting[sector] = 0
        f_ss_other_electricity[sector] = 0
        f_ss_other_gas[sector] = 0

        for region_name in data['lu_reg']:

            # HDD
            reg_hdd = ss_hdd_individ_region[region_name]
            reg_cdd = ss_cdd_individ_region[region_name]

            # Floor Area of sector
            reg_floor_area = data['ss_sector_floor_area_by'][region_name][sector]

            # Population
            reg_pop = data['population'][data['sim_param']['base_yr']][region_name]

            # National disaggregation factors
            f_ss_catering[sector] += reg_pop
            f_ss_computing[sector] += reg_pop
            f_ss_cooling_ventilation[sector] += reg_pop * reg_cdd
            f_ss_water_heating[sector] += reg_pop
            f_ss_space_heating[sector] += reg_floor_area * reg_hdd
            f_ss_lighting[sector] += reg_floor_area
            f_ss_other_electricity[sector] += reg_pop
            f_ss_other_gas[sector] += reg_pop

    # ---------------------------------------
    # Disaggregate according to enduse
    # ---------------------------------------
    for region_name in data['lu_reg']:
        ss_fueldata_disagg[region_name] = {}
        for sector in data['ss_sectors']:
            ss_fueldata_disagg[region_name][sector] = {}
            for enduse in data['ss_all_enduses']:

                # HDD
                reg_hdd = ss_hdd_individ_region[region_name]
                reg_cdd = ss_cdd_individ_region[region_name]

                # Floor Area of sector
                reg_floor_area = data['ss_sector_floor_area_by'][region_name][sector]

                # Population
                reg_pop = data['population'][data['sim_param']['base_yr']][region_name]

                if enduse == 'ss_catering':
                    reg_diasg_factor = reg_pop / f_ss_catering[sector]
                elif enduse == 'ss_computing':
                    reg_diasg_factor = reg_pop / f_ss_computing[sector]
                elif enduse == 'ss_cooling_ventilation':
                    reg_diasg_factor = (reg_pop * reg_cdd) / f_ss_cooling_ventilation[sector]
                elif enduse == 'ss_water_heating':
                    reg_diasg_factor = reg_pop / f_ss_water_heating[sector]
                elif enduse == 'ss_space_heating':
                    reg_diasg_factor = (reg_floor_area * reg_hdd) / f_ss_space_heating[sector]
                elif enduse == 'ss_lighting':
                    reg_diasg_factor = reg_floor_area / f_ss_lighting[sector]
                if enduse == 'ss_other_electricity':
                    reg_diasg_factor = reg_pop / f_ss_other_electricity[sector]
                elif enduse == 'ss_other_gas':
                    reg_diasg_factor = reg_pop / f_ss_other_gas[sector]

                # Disaggregate (fuel * factor)
                ss_fueldata_disagg[region_name][sector][enduse] = raw_fuel_sectors_enduses[sector][enduse] * reg_diasg_factor

    # TESTING Check if total fuel is the same before and after aggregation
    control_sum1 = 0
    control_sum2 = 0
    for reg in ss_fueldata_disagg:
        for sector in ss_fueldata_disagg[reg]:
            for enduse in ss_fueldata_disagg[reg][sector]:
                control_sum1 += np.sum(ss_fueldata_disagg[reg][sector][enduse])

    for sector in data['ss_sectors']:
        for enduse in data['ss_all_enduses']:
            control_sum2 += np.sum(raw_fuel_sectors_enduses[sector][enduse])

    #The loaded floor area must correspond to provided fuel sectors numers
    np.testing.assert_almost_equal(control_sum1, control_sum2, decimal=2, err_msg="")

    return ss_fueldata_disagg

def scrap_ts_disaggregate(data, fuel_national):
    """Disaggregate transport sector
    """
    fueldata_disagg = {}

    national_floorarea_sector = 0
    for region_name in data['lu_reg']:
        national_floorarea_sector += sum(data['ss_sector_floor_area_by'][region_name].values())

    # Iterate regions
    for region_name in data['lu_reg']:
        fueldata_disagg[region_name] = {}
        reg_floorarea_sector = sum(data['ss_sector_floor_area_by'][region_name].values())
        reg_disaggregation_factor = (1 / national_floorarea_sector) * reg_floorarea_sector

        reg_fuel_sector_enduse = reg_disaggregation_factor * fuel_national

        fueldata_disagg[region_name] = reg_fuel_sector_enduse

    return fueldata_disagg

def is_disaggregate(data, raw_fuel_sectors_enduses):
    """TODO: Disaggregate fuel for sector and enduses with floor
    area and GVA for sectors and enduses (IMPROVE)

    #TODO: DISAGGREGATE WITH OTHER DATA
    """
    is_fueldata_disagg = {}

    national_floorarea_sector = 0
    for region_name in data['lu_reg']:
        national_floorarea_sector += sum(data['ss_sector_floor_area_by'][region_name].values())

    # Iterate regions
    for region_name in data['lu_reg']:
        is_fueldata_disagg[region_name] = {}

        # Iterate sector
        for sector in data['is_sectors']:
            is_fueldata_disagg[region_name][sector] = {}

            # Sector specifid info
            reg_floorarea_sector = sum(data['ss_sector_floor_area_by'][region_name].values())

            # Iterate enduse
            for enduse in data['is_all_enduses']:
                national_fuel_sector_by = raw_fuel_sectors_enduses[sector][enduse]

                #print("national_fuel_sector_by: " + str(national_fuel_sector_by))
                # ----------------------
                # Disaggregating factors
                # TODO: IMPROVE. SHOW HOW IS DISAGGREGATED
                reg_disaggregation_factor = (1 / national_floorarea_sector) * reg_floorarea_sector

                # Disaggregated national fuel
                reg_fuel_sector_enduse = reg_disaggregation_factor * national_fuel_sector_by

                is_fueldata_disagg[region_name][sector][enduse] = reg_fuel_sector_enduse

    return is_fueldata_disagg

def rs_disaggregate(data, rs_national_fuel):
    """Disaggregate residential fuel demand

    Parameters
    ----------
    data : dict
        Data container
    rs_national_fuel : dict
        Fuel per enduse for residential submodel

    Returns
    -------
    rs_fueldata_disagg : dict
        Disaggregated fuel per enduse for every region (fuel[region][enduse])

    Note
    -----
    Used disaggregation factors for residential according
    to enduse (see Section XY Documentation TODO)
    """
    print("... disagreggate residential demand")

    # ---------------------------------------
    # Calculate heating degree days for regions
    # ---------------------------------------
    rs_hdd_individ_region = hdd_cdd.get_hdd_country(data['lu_reg'], data, 'rs_t_base_heating')

    total_pop = sum(data['population'][data['sim_param']['base_yr']].values()) # Total population
    # ---------------------------------------
    # Overall disaggregation factors per enduse
    # ---------------------------------------
    f_rs_lighting = 0
    fs_rs_cold = 0
    fs_rs_wet = 0
    fs_rs_consumer_electronics = 0
    fs_rs_home_computing = 0
    fs_rs_cooking = 0
    f_rs_space_heating = 0
    fs_rs_water_heating = 0

    for region_name in data['lu_reg']:

        # HDD
        reg_hdd = rs_hdd_individ_region[region_name]

        # Floor Area across all sectors
        reg_floor_area = data['rs_floorarea'][data['sim_param']['base_yr']][region_name]

        # Population
        reg_pop = data['population'][data['sim_param']['base_yr']][region_name]

        # National dissagregation factors
        f_rs_lighting += reg_floor_area
        fs_rs_cold += reg_pop
        fs_rs_wet += reg_pop
        fs_rs_consumer_electronics += reg_pop
        fs_rs_home_computing += reg_pop
        fs_rs_cooking += reg_pop
        f_rs_space_heating += reg_hdd * reg_floor_area
        fs_rs_water_heating += reg_pop

    # ---------------------------------------
    # Disaggregate according to enduse
    # ---------------------------------------
    rs_fueldata_disagg = {}
    for region_name in data['lu_reg']:
        rs_fueldata_disagg[region_name] = {}

        reg_pop = data['population'][data['sim_param']['base_yr']][region_name]
        reg_hdd = rs_hdd_individ_region[region_name]
        reg_floor_area = data['rs_floorarea'][data['sim_param']['base_yr']][region_name]

        # Disaggregate fuel depending on end_use
        for enduse in rs_national_fuel:
            if enduse == 'rs_lighting':
                reg_diasg_factor = reg_floor_area / f_rs_lighting
            elif enduse == 'rs_cold':
                reg_diasg_factor = reg_pop / fs_rs_cold
            elif enduse == 'rs_wet':
                reg_diasg_factor = reg_pop / fs_rs_wet
            elif enduse == 'rs_consumer_electronics':
                reg_diasg_factor = reg_pop / fs_rs_consumer_electronics
            elif enduse == 'rs_home_computing':
                reg_diasg_factor = reg_pop / fs_rs_home_computing
            elif enduse == 'rs_cooking':
                reg_diasg_factor = reg_pop / fs_rs_cooking
            elif enduse == 're_space_heating':
                reg_diasg_factor = (reg_hdd * reg_floor_area) / f_rs_space_heating
            elif enduse == 'rs_water_heating':
                reg_diasg_factor = reg_pop / fs_rs_water_heating
            else:
                reg_diasg_factor = reg_pop / total_pop

            # Disaggregate
            rs_fueldata_disagg[region_name][enduse] = rs_national_fuel[enduse] * reg_diasg_factor

    return rs_fueldata_disagg

def write_disagg_fuel(path_to_txt, data):
    """Write out disaggregated fuel

    Parameters
    ----------
    path_to_txt : str
        Path to txt file
    data : dict
        Data to write out
    """
    file = open(path_to_txt, "w")
    file.write("{}, {}, {}, {}".format(
        'region', 'enduse', 'fueltype', 'fuel') + '\n'
              )

    for region, enduses in data.items():
        for enduse, fuels in enduses.items():
            for fueltype, fuel in enumerate(fuels):
                file.write("{}, {}, {}, {}".format(
                    str.strip(region),
                    str.strip(enduse),
                    str(int(fueltype)),
                    str(float(fuel)) + '\n')
                          )
    file.close()

    return

def write_disagg_fuel_ts(path_to_txt, data):
    """Write out disaggregated fuel

    Parameters
    ----------
    path_to_txt : str
        Path to txt file
    data : dict
        Data to write out
    """
    file = open(path_to_txt, "w")
    file.write("{}, {}, {}".format(
        'region', 'fueltype', 'fuel') + '\n'
              )

    for region, fuels in data.items():
        for fueltype, fuel in enumerate(fuels):
            file.write("{}, {}, {}".format(
                str.strip(region), str(int(fueltype)), str(float(fuel)) + '\n')
                      )
    file.close()

    return

def write_disagg_fuel_sector(path_to_txt, data):
    """Write out disaggregated fuel

    Parameters
    ----------
    path_to_txt : str
        Path to txt file
    data : dict
        Data to write out
    """
    file = open(path_to_txt, "w")
    file.write("{}, {}, {}, {}, {}".format(
        'region', 'enduse', 'sector', 'fueltype', 'fuel') + '\n'
              )

    for region, sectors in data.items():
        for sector, enduses in sectors.items():
            for enduse, fuels in enduses.items():
                for fueltype, fuel in enumerate(fuels):
                    file.write("{}, {}, {}, {}, {}".format(
                        str.strip(region),
                        str.strip(enduse),
                        str.strip(sector),
                        str(int(fueltype)),
                        str(float(fuel)) + '\n')
                              )
    file.close()

    return

def run(path_main, local_data_path):
    """Function run script
    """
    print("... start script {}".format(os.path.basename(__file__)))

    # Load data and assumptions
    base_data = data_loader.load_paths(path_main, local_data_path)
    base_data = data_loader.load_fuels(base_data)
    base_data['assumptions'] = assumptions.load_assumptions(base_data)
    base_data['weather_stations'], temperature_data = data_loader.load_data_temperatures(
        os.path.join(base_data['paths']['path_scripts_data'], 'weather_data')
        )
    base_data['temperature_data'] = {}
    for weather_station, base_yr_temp in temperature_data.items():
        base_data['temperature_data'][weather_station] = {2015: base_yr_temp}

    # IMPROVE TODO: LOAD FLOOR AREA DATA
    # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
    base_data = data_loader.dummy_data_generation(base_data)
    # <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<

    # Disaggregation
    base_data = disaggregate_base_demand(base_data)

    #Write to csv file disaggregated demand
    write_disagg_fuel(
        os.path.join(path_main, 'data', 'data_scripts', 'disaggregated', 'rs_fueldata_disagg.csv'),
        base_data['rs_fueldata_disagg'])
    write_disagg_fuel_sector(
        os.path.join(path_main, 'data', 'data_scripts', 'disaggregated', 'ss_fueldata_disagg.csv'),
        base_data['ss_fueldata_disagg'])
    write_disagg_fuel_sector(
        os.path.join(path_main, 'data', 'data_scripts', 'disaggregated', 'is_fueldata_disagg.csv'),
        base_data['is_fueldata_disagg'])
    write_disagg_fuel_ts(
        os.path.join(path_main, 'data', 'data_scripts', 'disaggregated', 'ts_fueldata_disagg.csv'),
        base_data['ts_fueldata_disagg'])

    print("... finished script {}".format(os.path.basename(__file__)))
    return


text:
#!/usr/bin/env python3
# -*- coding:utf-8 -*-
# import the necessary package
import numpy as np
import cv2

# initialize the class labels and set the seed of the pseudorandom
# number generator so we can reproduce our results
labels = ["dog","cat"]
np.random.seed(1)

# randomly initialize our weight matrix and bias vector -- in a
# *real* training and classification task, these parameters would
# be *learned* by our model, but for the sake of this example
# let's use random values
W = np.random.randn(2,3072)
b = np.random.randn(2)

# load our example image, resize it, and then flatten it into our
# "feature vector" representation
orig = cv2.imread("../../../datasets/拉布拉多.jpg")
image = cv2.resize(orig, (32,32)).flatten()

# compute the output scores by taking the dot product between the
# weight matrix and image pixels, followed by adding in the bias
scores = W.dot(image) + b

# loop over the scores + labels and display them
for (label, score) in zip(labels, scores):
    print("[INFO] {}: {:.2f}".format(label, score))

# draw the label with the highest score on the image as our
# prediction
cv2.putText(orig, "Label: {}".format(labels[np.argmax(scores)]),
                                     (10,30), cv2.FONT_HERSHEY_SIMPLEX, 0.9, (0,255,0),2)

# display our input image
cv2.imshow("Image", orig)
cv2.waitKey(0)

max_stars_repo_name: ChyiYaqing/chyidlTutorial

max_issues_repo_name: ChyiYaqing/chyidlTutorial

max_forks_repo_name: ChyiYaqing/chyidlTutorial


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

df=pd.read_csv("corrected_data.csv")

# to get a new data frame with the desired columns
df_test=df[['drive-wheels','body-style','price']]

#to chane datatype of price column
df_test['price']=pd.to_numeric(df_test['price'],errors='coerce')

#to group the rows with drive-wheels and body-style and get the average price
df_grp = df_test.groupby(['drive-wheels','body-style'],as_index=False).mean()

#convert the data into a pivot table
df_pivot=df_grp.pivot(index='drive-wheels',columns='body-style')

#to make a heat map of the data
plt.pcolor(df_pivot,cmap='RdBu')
plt.colorbar()
plt.show()

max_stars_repo_name: Arya-A-Nair/learning-data-science

max_issues_repo_name: Arya-A-Nair/learning-data-science

max_forks_repo_name: Arya-A-Nair/learning-data-science


text:
import numpy as np
from tqdm import tqdm

def ssnmf_core(V, rdim, iter_num, W0, H0):
    if (V < 0).any():
        raise ValueError("Negative values in data!")

    vdim = V.shape[0]
    samples = V.shape[1]

    W = W0
    H = H0

    for iteration in range(1, iter_num + 1):

        # Update using standard NMF multiplicative update rule
        H = H * (W.T @ V) / (W.T @ W @ H + 1e-9)

        # Re-normalize so that rows of H have constant energy
        norms = np.sqrt(np.sum(H.T ** 2, 0))
        H = H / np.outer(norms.T, np.ones((1, samples)))
        W = W * np.outer(np.ones((vdim, 1)), norms)

        # Update using standard NMF multiplicative update rule
        W = W * (V @ H.T) / (W @ H @ H.T + 1e-9)

    return (W, H)

def ssnmf(spectrogramData2D, number_basis, number_iter):
    '''
    2022-01-30 (v01)
        - Fuh-Cherng (Fuh) Jeng converted part of the nmfsc(), that was originally written by Tzu-Hao (Harry) Lin in Matlab, to Python
        - Fuh made a few minor changes
    2022-02-22 (v02)
        - Fuh modified this script so that the input spectrogramData2D must be a 2D matrix. Otherwise, this script will not work.
        - As such, no reshape() is needed within this function.
        - In other words, any reshape() must be performed prior to feeding spectrogramData2D into this function.
    '''

    data = 20 * np.log10(spectrogramData2D)   # convert nV to dB (adopted from Lin et al., 2017).
    # IMPORTANT: This converion is just for the purpose of optimizing W and H. It has no effect on the amplitude unit of ssnmf() output. The unit of ssnmf() output will always be the same the unit of the input data. For this study, the unit of the input data is nV.

    data[data < 0] = 0

    n = data.shape[1]
    weight = np.ones((1, n))

    W = abs(np.random.randn(data.shape[0], number_basis))
    H = abs(np.random.randn(number_basis, n))
    H = H / np.outer(np.sqrt(np.sum(H ** 2, 1)), np.ones((1, n)))

    # NMF decomposition
    print('performing ssnmf(), please wait...', flush=True)
    for run in tqdm(range(number_iter)):
        H[0, :] = weight * np.sum(H[0, :]) / np.sum(weight)
        (W, H) = ssnmf_core(data, number_basis, 1, W, H)

    output = np.zeros((data.shape[0], data.shape[1], number_basis))
    np.seterr(invalid='ignore')
    for run in range(number_basis):
        output[:, :, run] = spectrogramData2D * np.outer(W[:, run], H[run, :]) / (W @ H)

    np.seterr(invalid='warn')
    output[np.isnan(output)] = 0

    return (output, W, H)


text:
import numpy as np
from scipy.stats import truncnorm

from prml.nn.array.array import asarray

def normal(mean, std, size):
    return asarray(np.random.normal(mean, std, size))

def truncnormal(min, max, scale, size):
    return asarray(truncnorm(a=min, b=max, scale=scale).rvs(size))

max_stars_repo_name: alexandru-dinu/PRML

max_issues_repo_name: alexandru-dinu/PRML

max_forks_repo_name: alexandru-dinu/PRML

reason: import numpy,from scipy


text:
import argparse
import os
import time
import matplotlib.pyplot as plt
import numpy as np
import torch
import torch.utils.data
import torchvision.utils as vutils
import pandas as pd
from math import pi, sin, cos, sqrt, log
from toolz import partial, curry
from torch import Tensor
from torch import nn, optim, distributions
from torchvision import datasets, transforms, models
from torchvision.utils import save_image, make_grid
from typing import Callable, Iterator, Union, Optional, TypeVar
from typing import List, Set, Dict, Tuple
from typing import Mapping, MutableMapping, Sequence, Iterable
from typing import Union, Any, cast
# my own sauce
from my_torch_utils import denorm, normalize, mixedGaussianCircular
from my_torch_utils import fclayer, init_weights
from my_torch_utils import plot_images, save_reconstructs, save_random_reconstructs
from my_torch_utils import fnorm, replicate, logNorm
from my_torch_utils import scsimDataset

import scsim.scsim as scsim

device = "cuda" if torch.cuda.is_available() else "cpu"
dataSet = scsimDataset("data/scrnasim/counts.npz",
        "data/scrnasim/cellparams.npz")
trainD, testD = dataSet.__train_test_split__(8500)

trainLoader = torch.utils.data.DataLoader(trainD, batch_size=128, shuffle=True)
testLoader = torch.utils.data.DataLoader(testD, batch_size=128, shuffle=True)

xs, ls = trainLoader.__iter__().next()

class Classifier(nn.Module):
    def __init__(self, nin=10**4, nh=10**3, nclasses=10):
        super(Classifier, self).__init__()
        self.nclasses = nclasses
        self.nin = nin
        self.main = nn.Sequential(
                fclayer(nin=nin, nout=nh, batchnorm=True, dropout=0.2,
                    activation=nn.LeakyReLU(),),
                nn.Linear(nh, nclasses),
                nn.LogSoftmax(dim=1),
                )
        return

    def forward(self,x):
        logp = self.main(x)
        return logp

model = Classifier()

model(xs).shape

ys = model(xs)

criterion = nn.NLLLoss(reduction="mean")

criterion(ys, ls-1)

model.to(device)

model.apply(init_weights)

optimizer = optim.Adam(model.parameters())

for epoch in range(15):
    for idx, (data, labels) in enumerate(trainLoader):
        x = data.to(device)
        target = (labels - 1).to(device)
        model.train()
        model.zero_grad()
        logprob = model(x)
        loss = criterion(logprob, target)
        loss.backward()
        optimizer.step()
        if idx % 100 == 0:
            print(
                "loss:\n",
                loss.item(),
                )

# testing
xs, ls = testLoader.__iter__().next()

# convert the labels to 0-indexing
ls -= 1

model.to("cpu")

ys = model(xs)

probs = ys.exp()

criterion(ys, ls)

predicts = probs.argmax(axis=1)

predicts
ls

### generate new simulated dataset
simulator = scsim.scsim(ngenes=5000, nproggenes=300, ncells=10000,
        ngroups=10, seed=42,libloc=7.64, libscale=0.78,
        mean_rate=7.68, mean_shape=0.34, expoutprob=0.00286, expoutloc=6.15,
        expoutscale=0.49,
        )

ngenes = 10**4
ncells = 10**4
K=10
deprob = .025
progdeloc = deloc = deval = 1
descale = 1.0
progcellfrac = .35
deprob = .025
doubletfrac = .0
ndoublets=int(doubletfrac*ncells)
nproggenes = 400
nproggroups = int(K/3)
proggroups = list(range(1, nproggroups+1))
randseed = 42

simulator = scsim.scsim(
    ngenes=ngenes,
    ncells=ncells,
    ngroups=K,
    libloc=7.64,
    libscale=0.78,
    mean_rate=7.68,
    mean_shape=0.34,
    expoutprob=0.00286,
    expoutloc=6.15,
    expoutscale=0.49,
    diffexpprob=deprob,
    diffexpdownprob=0.0,
    diffexploc=deloc,
    diffexpscale=descale,
    bcv_dispersion=0.448,
    bcv_dof=22.087,
    ndoublets=ndoublets,
    nproggenes=nproggenes,
    progdownprob=0.0,
    progdeloc=progdeloc,
    progdescale=descale,
    progcellfrac=progcellfrac,
    proggoups=proggroups,
    minprogusage=0.1,
    maxprogusage=0.7,
    seed=randseed,

)

simulator.simulate()

simulator.cellparams.to_csv(path_or_buf="data/scrnasim/cellparams2.tsv.gz",
        sep='\t',
        compression="gzip",
        )

df = pd.read_csv("./data/scrnasim/cellparams2.tsv.gz",
        sep='\t', index_col=0, compression='gzip',
        )

df

f = np.load("data/scrnasim/cellparams2.tsv", allow_pickle=True)

df = scsim.load_df("data/scrnasim/cellparams.npz")

df

scsim.save_df(df, "data/scrnasim/blabla.gz")

df = scsim.load_df("data/scrnasim/blabla.gz.npz")

df

max_stars_repo_name: zelhar/MPGVAE

max_issues_repo_name: zelhar/MPGVAE

max_forks_repo_name: zelhar/MPGVAE