text:
import os
import csv

from opencmiss.utils.zinc.general import ChangeManager
from opencmiss.zinc.context import Context
from opencmiss.zinc.field import Field, FieldFindMeshLocation
from opencmiss.zinc.node import Node
from opencmiss.utils.zinc.finiteelement import getMaximumNodeIdentifier
from opencmiss.zinc.element import Element, Elementbasis
from opencmiss.zinc.result import RESULT_OK
from opencmiss.utils.zinc.finiteelement import findNodeWithName
from opencmiss.utils.zinc.field import findOrCreateFieldGroup, findOrCreateFieldNodeGroup, findOrCreateFieldCoordinates
from opencmiss.utils.zinc.field import findOrCreateFieldStoredMeshLocation
from transformation import read_input_target, align_mesh, convert_to_transformation_matrix, write_to_exnode
import numpy as np

class CsvToExnode:
    def __init__(self, filepath, output, nodeOffset, group_name, coordinates_name='coordinates'):
        self._fileName = os.path.basename(filepath)
        self._dirname = os.path.dirname(filepath)
        self._output = output
        self._inputFile = filepath
        self._nodeOffset = nodeOffset
        self._coordinates_name = coordinates_name
        self._group_name = group_name
        self._context = Context("CsvToExnode")
        self._region = self._context.getDefaultRegion()
        self._fieldmodule = self._region.getFieldmodule()
        self._onlyCoordinates = True

    def csv_to_exnode(self):
        self.convert_csv_to_exf()
        self.convert_exf_to_exnode()

    def convert_csv_to_exf(self):
        coordinates = findOrCreateFieldCoordinates(self._fieldmodule, components_count=3)
        cache = self._fieldmodule.createFieldcache()

        #################
        # Create nodes
        #################

        nodes = self._fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES)
        nodetemplate = nodes.createNodetemplate()
        nodetemplate.defineField(coordinates)
        nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_VALUE, 1)
        if not self._onlyCoordinates:
            nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS1, 1)

        nodeIdentifier = self._nodeOffset
        with open(self._inputFile, 'r') as csvf:
            reader = csv.reader(csvf)
            next(reader, None)
            for point in reader:
                if nodeIdentifier == self._nodeOffset:
                    x0 = [float(p) for p in point[0:3]]
                    nodeIdentifier += 1
                elif nodeIdentifier == self._nodeOffset + 1:
                    x1 = [float(p) for p in point[0:3]]
                    node = nodes.createNode(nodeIdentifier - 1, nodetemplate)
                    cache.setNode(node)
                    coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, x0)
                    if not self._onlyCoordinates:
                        dx_ds1 = [x1[c] - x0[c] for c in range(3)]
                        coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, dx_ds1)
                    node = nodes.createNode(nodeIdentifier, nodetemplate)
                    cache.setNode(node)
                    coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, x1)
                    if not self._onlyCoordinates:
                        dx_ds1 = [x1[c] - x0[c] for c in range(3)]
                        coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, dx_ds1)
                    x0 = [d for d in x1]
                    nodeIdentifier = nodeIdentifier + 1
                else:
                    x1 = [float(p) for p in point[0:3]]
                    node = nodes.createNode(nodeIdentifier, nodetemplate)
                    cache.setNode(node)
                    coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, x1)
                    if not self._onlyCoordinates:
                        dx_ds1 = [x1[c] - x0[c] for c in range(3)]
                        coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, dx_ds1)
                    x0 = [d for d in x1]
                    nodeIdentifier = nodeIdentifier + 1

        outputFile = os.path.join(self._dirname, self._fileName.split('.')[0]+'.exf')
        self._region.writeFile(outputFile)

    def convert_exf_to_exnode(self, **kwargs):
        """
        Write model nodes to file.
        """
        filepath = kwargs['filepath'] if 'filepath' in kwargs else os.path.join(self._dirname, self._fileName.split('.')[0] + '.exf')
        filename = os.path.basename(filepath)
        context = Context('convert to exnode')
        region = context.getDefaultRegion()
        region.readFile(filepath)
        fieldmodule = region.getFieldmodule()
        coordinates = fieldmodule.findFieldByName('coordinates').castFiniteElement()
        coordinates.setName(self._coordinates_name)
        nodes = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES)

        markerGroup = findOrCreateFieldGroup(fieldmodule, self._group_name)
        markerPoints = findOrCreateFieldNodeGroup(markerGroup, nodes).getNodesetGroup()
        is_all = fieldmodule.createFieldConstant(1)
        markerPoints.addNodesConditional(is_all)

        sir = region.createStreaminformationRegion()
        sir.setRecursionMode(sir.RECURSION_MODE_OFF)
        srm = sir.createStreamresourceMemory()
        sir.setResourceGroupName(srm, self._group_name)
        sir.setResourceDomainTypes(srm, Field.DOMAIN_TYPE_NODES)
        result = region.write(sir)
        assert result == RESULT_OK
        exnodeFileName = filename.split('.')[0] + ".exnode"
        working_directory = os.path.dirname(filepath)
        exnodeFilepath = os.path.join(working_directory, exnodeFileName)
        region.writeFile(exnodeFilepath)

class CsvToExf:
    def __init__(self, file, nodeOffset):
        fileName = os.path.basename(file)
        dirname = os.path.dirname(file)
        self._inputFile = file
        self._context = Context("CsvToEx")
        self._region = self._context.getDefaultRegion()
        self._fieldmodule = self._region.getFieldmodule()
        self._onlyCoordinates = True

        coordinates = findOrCreateFieldCoordinates(self._fieldmodule, components_count=3)
        cache = self._fieldmodule.createFieldcache()

        #################
        # Create nodes
        #################

        nodes = self._fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES)
        nodetemplate = nodes.createNodetemplate()
        nodetemplate.defineField(coordinates)
        nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_VALUE, 1)
        if not self._onlyCoordinates:
            nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS1, 1)

        nodeIdentifier = nodeOffset
        with open(self._inputFile, 'r') as csvf:
            reader = csv.reader(csvf)
            next(reader, None)
            for point in reader:
                if nodeIdentifier == nodeOffset:
                    x0 = [float(p) for p in point[0:3]]
                    nodeIdentifier += 1
                elif nodeIdentifier == nodeOffset + 1:
                    x1 = [float(p) for p in point[0:3]]
                    node = nodes.createNode(nodeIdentifier - 1, nodetemplate)
                    cache.setNode(node)
                    coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, x0)
                    if not self._onlyCoordinates:
                        dx_ds1 = [x1[c] - x0[c] for c in range(3)]
                        coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, dx_ds1)
                    node = nodes.createNode(nodeIdentifier, nodetemplate)
                    cache.setNode(node)
                    coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, x1)
                    if not self._onlyCoordinates:
                        dx_ds1 = [x1[c] - x0[c] for c in range(3)]
                        coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, dx_ds1)
                    x0 = [d for d in x1]
                    nodeIdentifier = nodeIdentifier + 1
                else:
                    x1 = [float(p) for p in point[0:3]]
                    node = nodes.createNode(nodeIdentifier, nodetemplate)
                    cache.setNode(node)
                    coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, x1)
                    if not self._onlyCoordinates:
                        dx_ds1 = [x1[c] - x0[c] for c in range(3)]
                        coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, dx_ds1)
                    x0 = [d for d in x1]
                    nodeIdentifier = nodeIdentifier + 1

        outputFile = os.path.join(dirname, fileName.split('.')[0]+'.exf')
        self._region.writeFile(outputFile)

class EmbedOrgan:
    """
    A class for extracting organ marker coordinates and embed the fiducial markers in the body-coordinates
    """

    def __init__(self, organ_name, organ_filepath, working_directory, *body_filepath):
        self._organ_name = organ_name
        self._working_directory = working_directory
        self._organ_filepath = organ_filepath
        self._body_filepath = body_filepath
        self._body_exnodefilepath = body_filepath[0]
        self._organ_fiducial_markers = None

        # self.get_organ_fiducials(organ_filename)
        # context = Context('extract_marker_coordinates')
        # self._region = context.getDefaultRegion()
        # self._organ_dirname = organ_dirname
        # organ_file = organ_filename
        # self._region.readFile(organ_file)
        # self._fieldmodule = self._region.getFieldmodule()
        # self._mesh = self.get_highest_dimension_mesh()
        # self._nodes = self._fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES)
        # self._coordinates = self._fieldmodule.findFieldByName('coordinates').castFiniteElement()
        # self._markerNodes = None

        # organ_reference.exf
        # markerGroup = self._fieldmodule.findFieldByName("marker")
        # if markerGroup.isValid():
        #     markerGroup = markerGroup.castGroup()
        #     markerNodeGroup = markerGroup.getFieldNodeGroup(self._nodes)
        #     if markerNodeGroup.isValid():
        #         self._markerNodes = markerNodeGroup.getNodesetGroup()
        #
        # if self._markerNodes and (self._markerNodes.getSize() > 0):
        #     filename = os.path.basename(organ_filename)
        #     markerFilename = filename.split('.')[0] + "_marker.csv"
        #     markerFilePath = os.path.join(organ_dirname, markerFilename)
        #     # markerFilename = os.path.splitext(organ_file)[0] + "_marker.csv"
        #     with open(markerFilePath, 'w') as outstream:
        #         self._writeMarkers(outstream)

        # self.read_landmarks()
        # self.get_marker_location()
        # markerGroup2 = self._fieldmodule.findFieldByName("marker").castGroup()
        # markerNodes = markerGroup2.getFieldNodeGroup(self._nodes).getNodesetGroup()
        # markerName = self._fieldmodule.findFieldByName("marker_name")
        # assert markerName.isValid()
        # markerLocation = self._fieldmodule.findFieldByName("marker_location")
        # markerLocation2 = self._fieldmodule.createFieldFindMeshLocation(self._markerCoordinates, self._coordinates, self._mesh)
        # a = markerLocation2.setManaged(True)
        # markerLocation2.setName("Elias")
        # # assert self._markerLocation.isValid(), 'invlaide marker location'
        # # test apex marker point
        # cache = self._fieldmodule.creorgan_filepathateFieldcache()
        # node = findNodeWithName(markerNodes, markerName, "apex of heart")
        # assert node.isValid(), 'node is invalid'
        # cache.setNode(node)
        # nodeId = node.getIdentifier()
        # element, xi = markerLocation2.evaluateMeshLocation(cache, 3)
        # elementId = element.getIdentifier()
        # self._region.writeFile('rat_body_data500.exf')
        a=1

    def get_maximum_node_identifier(self):
        context = Context('node number')
        region = context.getDefaultRegion()
        region.readFile(self._body_exnodefilepath)
        fieldmodule = region.getFieldmodule()
        coordinates = fieldmodule.findFieldByName('coordinates').castFiniteElement()
        nodes = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES)
        self._maximum_nodeId = getMaximumNodeIdentifier(nodes)
        return self._maximum_nodeId

    def set_organ_filepath(self, organ_filepath):
        self._organ_filepath = organ_filepath

    def get_marker_coordinates(self, outputpath, *filepath, **kwargs):
        working_directory = kwargs['working_directory'] if 'working_directory' in kwargs else self._working_directory
        if 'marker_list' in kwargs:
            marker_list = kwargs['marker_list']
        else:
            marker_list = []

        context = Context('organ_fiducial_coordinates')
        region = context.getDefaultRegion()
        for f in filepath:
            region.readFile(f)
        fieldmodule = region.getFieldmodule()
        coordinates = fieldmodule.findFieldByName('coordinates').castFiniteElement()
        nodes = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES)

        coordinatesCount = coordinates.getNumberOfComponents()
        cache = fieldmodule.createFieldcache()

        markerLocation = fieldmodule.findFieldByName("marker_location")
        markerName = fieldmodule.findFieldByName("marker_name")
        markerGroup = fieldmodule.findFieldByName("marker")

        markers = []
        if markerGroup.isValid():
            markerGroup = markerGroup.castGroup()
            markerNodeGroup = markerGroup.getFieldNodeGroup(nodes)
            if markerNodeGroup.isValid():
                markerNodes = markerNodeGroup.getNodesetGroup()

        if markerLocation.isValid() and markerName.isValid():
            with ChangeManager(fieldmodule):
                markerCoordinates = fieldmodule.createFieldEmbedded(coordinates, markerLocation)
                nodeIter = markerNodes.createNodeiterator()
                node = nodeIter.next()
                with open(outputpath, 'w') as outstream:
                    outstream.write('x,y,z,landmark\n')
                    while node.isValid():
                        cache.setNode(node)
                        result, x = markerCoordinates.evaluateReal(cache, coordinatesCount)
                        if result == RESULT_OK:
                            name = markerName.evaluateString(cache)
                            if (name in marker_list) or (marker_list == []):
                                markers.append({'landmark': name, 'x': x})
                                outstream.write(",".join(str(s) for s in x) + "," + name + "\n")
                        node = nodeIter.next()

        return markers

    def get_organ_fiducials(self, organ_filename):
        context = Context('extract_marker_coordinates')

    def read_landmarks(self, file):
        # file = r'C:\Users\egha355\Desktop\embed\modify_heart\whole-heart2.exf'
        # with open(file, 'r') as lmf:
        #     reader = csv.reader(lmf)
        #     next(reader, None)
        #     landmark_name = []
        #     landmark_coordinates = []
        #     for point in reader:
        #         landmark_name.append(point[0])
        #         landmark_coordinates.append([float(point[i]) for i in range(1, 4)])
        #     self._landmark_name = landmark_name
        #     self._landmark_coordinates = landmark_coordinates

        with open(file, 'r') as fsc:
            reader = csv.reader(fsc)
            next(reader, None)
            markers = []
            for landmark in reader:
                markers.append({'landmark': landmark[-1], 'x':landmark[:-1]})

        self._organ_landmarks = markers

    def get_highest_dimension_mesh(self):
        for dimension in range(3, 1, -1):
            mesh = self._fieldmodule.findMeshByDimension(dimension)
            if mesh.getSize() > 0:
                break
        return mesh

    # def _writeMarkers(self, outstream):
    #     coordinatesCount = self._coordinates.getNumberOfComponents()
    #     cache = self._fieldmodule.createFieldcache()
    #     # cache2 = self._fieldmodule.createFieldcache()
    #     markerLocation = self._fieldmodule.findFieldByName("marker_location")
    #     markerName = self._fieldmodule.findFieldByName("marker_name")
    #
    #     # organ_marker_coordinates = self._fieldmodule.createFieldFiniteElement(3)
    #     # nodetemplate = self._markerNodes.createNodetemplate()
    #     # nodetemplate.defineField(organ_marker_coordinates)
    #     # nodetemplate.setValueNumberOfVersions(organ_marker_coordinates, -1, Node.VALUE_LABEL_VALUE, 1)
    #     # nodeIdentifier = 8000100
    #
    #     if markerLocation.isValid() and markerName.isValid():
    #         with ChangeManager(self._fieldmodule):
    #             markerCoordinates = self._fieldmodule.createFieldEmbedded(self._coordinates, markerLocation)
    #             nodeIter = self._markerNodes.createNodeiterator()
    #             node = nodeIter.next()
    #             while node.isValid():
    #                 cache.setNode(node)
    #                 result, x = markerCoordinates.evaluateReal(cache, coordinatesCount)
    #
    #                 # node2 = self._markerNodes.createNode(nodeIdentifier, nodetemplate)
    #                 # cache2.setNode(node2)
    #                 # organ_marker_coordinates.setNodeParameters(cache2, -1, Node.VALUE_LABEL_VALUE, 1, x)
    #                 if result == RESULT_OK:
    #                     name = markerName.evaluateString(cache)
    #                     outstream.write(",".join(str(s) for s in x) + "," + name + "\n")
    #                 node = nodeIter.next()
    #
    #                 # nodeIdentifier = nodeIdentifier + 1
    #             self._markerCoordinates = markerCoordinates

    # def _write_marker_locations(self, outstream):
    #     coordinatesCount = self._coordinates.getNumberOfComponents()
    #     cache = self._fieldmodule.createFieldcache()
    #     markerLocation = self._fieldmodule.findFieldByName("marker_location")
    #     markerName = self._fieldmodule.findFieldByName("marker_name")
    #     if markerLocation.isValid() and markerName.isValid():
    #         with ChangeManager(self._fieldmodule):
    #             markerCoordinates = self._fieldmodule.createFieldEmbedded(self._coordinates, markerLocation)
    #             nodeIter = self._markerNodes.createNodeiterator()
    #             node = nodeIter.next()
    #             while node.isValid():
    #                 cache.setNode(node)
    #                 result, x = markerCoordinates.evaluateReal(cache, coordinatesCount)
    #                 if result == RESULT_OK:
    #                     name = markerName.evaluateString(cache)
    #                     outstream.write(",".join(str(s) for s in x) + "," + name + "\n")
    #                 node = nodeIter.next()
    #             self._markerCoordinates = markerCoordinates.castFiniteElement()

    # def get_marker_location(self):
    #     self._markerLocation = self._fieldmodule.createFieldFindMeshLocation(self._markerCoordinates, self._coordinates, self._mesh)
    #     assert RESULT_OK == self._markerLocation.setSearchMesh(self._mesh)
    #     self._markerLocation.setSearchMode(FieldFindMeshLocation.SEARCH_MODE_EXACT)

    def get_element_xi(self, bodyNodeFile, bodyElemFile, organ_finalstate):
        context = Context('convert to exnode')
        region = context.getDefaultRegion()

        # for f in body_filepath:
        #     region.readFile(f)
        region.readFile(bodyNodeFile)
        region.readFile(bodyElemFile)
        region.readFile(organ_finalstate)

        fieldmodule = region.getFieldmodule()
        mesh3d = fieldmodule.findMeshByDimension(3)
        coordinates = fieldmodule.findFieldByName('coordinates').castFiniteElement()
        nodes = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES)
        # TODO make it possible to get the bodyFile instead of this. Note becasue organ file has markers the body file should not have any marker.
        organ_marker_coordinates = fieldmodule.findFieldByName(organ_name + '_marker_coordinates').castFiniteElement()
        organ_locations = fieldmodule.createFieldFindMeshLocation(organ_marker_coordinates, coordinates, mesh3d)
        organ_locations.setSearchMode(FieldFindMeshLocation.SEARCH_MODE_EXACT)

        # markerGroup = fieldmodule.findFieldByName("marker").castGroup()
        # markerNodes = markerGroup.getFieldNodeGroup(nodes).getNodesetGroup()
        # markerName = fieldmodule.findFieldByName("marker_name")

        organMarkerGroup = fieldmodule.findFieldByName(organ_name + '_marker').castGroup()
        organMarkerNodesetGroup = organMarkerGroup.getFieldNodeGroup(nodes).getNodesetGroup()
        # markerName = fieldmodule.findFieldByName("marker_name")

        cache = fieldmodule.createFieldcache()
        nodeIter = organMarkerNodesetGroup.createNodeiterator()
        node = nodeIter.next()

        element_xi_filepath = os.path.join(self._working_directory, os.path.basename(organ_finalstate).split('_')[0] + '_elementxi.csv')

        with open(element_xi_filepath, 'w', newline='') as elxif:
            writer = csv.writer(elxif)
            writer.writerow(['marker', 'element', 'xi1', 'xi2', 'xi3'])
            co = 0
            while node.isValid():
                cache.setNode(node)
                element, xi = organ_locations.evaluateMeshLocation(cache, 3)
                print('Node: ', str(co+403)+'\n', element.getIdentifier(), xi[0], xi[1], xi[2], '\n', '"'+self._organ_landmarks[co]['landmark']+'"')
                writer.writerow([self._organ_landmarks[co]['landmark'], element.getIdentifier(), xi[0], xi[1], xi[2]])
                co += 1
                node = nodeIter.next()

        # How to add the new marker locations?
        # markerGroup = fieldmodule.findFieldByName("marker_location").castGroup()
        # markerNodes = markerGroup.getFieldNodeGroup(nodes).getNodesetGroup()
        # markerName = fieldmodule.findFieldByName("marker_name")

        # get marker names

    # def writeReference(self, organ_filename, organ_name):
    #     context = Context('organ_reference')
    #     region = context.getDefaultRegion()
    #     fieldmodule = region.getFieldmodule()
    #     region.readFile(organ_filename)
    #     coordinates = fieldmodule.findFieldByName('coordinates').castFiniteElement()
    #     coordinatesCount = coordinates.getNumberOfComponents()
    #     cache = fieldmodule.createFieldcache()
    #     markerLocation = fieldmodule.findFieldByName("marker_location")
    #     markerName = fieldmodule.findFieldByName("marker_name")
    #
    #     nodes = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES)
    #
    #     # organ_reference.exf
    #     markerGroup = fieldmodule.findFieldByName("marker")
    #     if markerGroup.isValid():
    #         markerGroup = markerGroup.castGroup()
    #         markerNodeGroup = markerGroup.getFieldNodeGroup(nodes)
    #         if markerNodeGroup.isValid():
    #             markerNodes = markerNodeGroup.getNodesetGroup()
    #     # organ_name = 'heart'
    #     working_directory = os.path.join(os.path.dirname(organ_filename), organ_name)
    #     filepath = os.path.join(working_directory, organ_name + '_reference.csv')
    #
    #     with open(filepath, 'w', newline='') as reff:
    #         writer = csv.writer(reff)
    #         writer.writerow(['x', 'y', 'z', 'landmark'])
    #         if markerLocation.isValid() and markerName.isValid():
    #             with ChangeManager(fieldmodule):
    #                 markerCoordinates = fieldmodule.createFieldEmbedded(coordinates, markerLocation)
    #                 nodeIter = markerNodes.createNodeiterator()
    #                 node = nodeIter.next()
    #                 while node.isValid():
    #                     cache.setNode(node)
    #                     result, x = markerCoordinates.evaluateReal(cache, coordinatesCount)
    #
    #                     if result == RESULT_OK:
    #                         name = markerName.evaluateString(cache)
    #                         x.append(name)
    #                         writer.writerow(x)
    #                         # outstream.write(",".join(str(s) for s in x) + "," + name + "\n")
    #                     node = nodeIter.next()

    # def findCoordInOtherTimes():
    #     # First read the files with different times.
    #     context = Context('Find_coordinates')
    #     region = context.getDefaultRegion()
    #     timeKeeper = context.getTimekeepermodule()
    #     sir = region.createStreaminformationRegion()
    #     sir.setAttributeReal(sir.ATTRIBUTE_TIME, 0.0)
    #     result = region.read(sir)
    #     assert result == RESULT_OK
    #     parentDirectory = r'C:\Users\egha355\Desktop\embed\modify_heart'
    #     filename = 'bodycoordAligned_node4.exf'
    #     filepath = os.path.join(parentDirectory, filename)
    #     region.readFile(filepath)
    #     sir.setAttributeReal(sir.ATTRIBUTE_TIME, 1.0)
    #     a=1
    #
    # # Find heart marker coordinates in different time steps.
    # # Fieldmodule::createFieldTimeLookup	(	const Field & 	sourceField,
    # # const Field & 	timeField
    # # )
    # # :Fieldmodule::createFieldConstant	(	int 	valuesCount,
    # # const double * 	valuesIn
    # # )
    # # gfx def field marker_coordinates embedded field coordinates element_xi marker_location
    # # gfx define field heart/reftime constant 3.0;
    # # gfx define field heart/coordinates_reftime coordinate_system rectangular_cartesian time_lookup field coordinates time_field reftime;

    # def findCoordInOtherTimes2(self, exnode, exelem, organ_name, parentDirectory, working_directory, target_filename, final_statecsv):
    #     context = Context('Find_coordinates')
    #     region = context.getDefaultRegion()
    #     fieldmodule = region.getFieldmodule()
    #     # parentDirectory = r'C:\Users\egha355\Desktop\embed\modify_heart'
    #     markerFilename = 'rat1_fit1_node_static.exf'
    #     markerFilepath = os.path.join(parentDirectory, markerFilename)
    #
    #     result0 = region.readFile(exnode)
    #     result1 = region.readFile(exelem)
    #     result2 = region.readFile(markerFilepath)
    #
    #     coordinates = fieldmodule.findFieldByName('coordinates').castFiniteElement()
    #     coordinatesCount = coordinates.getNumberOfComponents()
    #     cache = fieldmodule.createFieldcache()
    #     markerLocation = fieldmodule.findFieldByName("marker_location")
    #     markerName = fieldmodule.findFieldByName("marker_name")
    #
    #     nodes = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES)
    #
    #     with open(final_statecsv, 'r') as fsc:
    #         reader = csv.reader(fsc)
    #         next(reader, None)
    #         marker_name = []
    #         for landmark in reader:
    #             marker_name.append(landmark[-1])
    #
    #     # organ_reference.exf
    #     markerGroup = fieldmodule.findFieldByName("marker")
    #     if markerGroup.isValid():
    #         markerGroup = markerGroup.castGroup()
    #         markerNodeGroup = markerGroup.getFieldNodeGroup(nodes)
    #         if markerNodeGroup.isValid():
    #             markerNodes = markerNodeGroup.getNodesetGroup()
    #     # working_directory = os.path.join(parentDirectory, organ_name)
    #     # filepath = os.path.join(working_directory, 'trans_{}.csv'.format(organ_name))
    #     filepath = target_filename
    #
    #     with open(filepath, 'w', newline='') as trf:
    #         writer = csv.writer(trf)
    #         writer.writerow(['x', 'y', 'z', 'landmark'])
    #         if markerLocation.isValid() and markerName.isValid():
    #             with ChangeManager(fieldmodule):
    #                 markerCoordinates = fieldmodule.createFieldEmbedded(coordinates, markerLocation)
    #                 nodeIter = markerNodes.createNodeiterator()
    #                 node = nodeIter.next()
    #                 while node.isValid():
    #                     cache.setNode(node)
    #                     result, x = markerCoordinates.evaluateReal(cache, coordinatesCount)
    #
    #                     if result == RESULT_OK:
    #                         name = markerName.evaluateString(cache)
    #                         if name in marker_name:
    #                             x.append(name)
    #                             writer.writerow(x)
    #                         # outstream.write(",".join(str(s) for s in x) + "," + name + "\n")
    #                     node = nodeIter.next()

    def compute_transformation_matrix(self, working_directory, reference_filename, target_filename):
        # working_directory = r'C:\Users\egha355\Desktop\embed\modify_heart\heart'
        # reference_filename = 'heart_reference.csv'
        # target_filename = 'trans_heart0.csv'
        reference_filepath = os.path.join(working_directory, reference_filename)
        target_filepath = os.path.join(working_directory, target_filename)
        node = 601

        output_file =os.path.join(working_directory, os.path.basename(target_filename).split('.')[0] + '.exnode')

        xinput, xtarget, mean_input = read_input_target(reference_filepath, target_filepath)
        numberOfPoints = len(xinput)
        d, Z, tform = align_mesh(xtarget, xinput, scaling=True, reflection='best')
        print(d, Z)
        print(tform)
        transform, translate = convert_to_transformation_matrix(tform)
        print('transform\n', transform)
        print('translate\n', translate)
        print('s*Y*R+c')
        print(np.matmul(np.array(xinput),tform['scale']*tform['rotation'])+tform['translation'])
        print('xtarget\n',xtarget)
        print('xinputplus1 * (transform+translate.T)')
        print(np.matmul(np.concatenate((xinput, np.ones((numberOfPoints,1))), axis=1), transform+translate.T))
        write_to_exnode(node, transform, translate, mean_input, output_file, tform)

    def embed_organ(self, organ_finalstatepath, organ_filepath, bodyNodeFile, bodyElemFile):
        nodeOffset = self.get_maximum_node_identifier() + 100
        organ_finalstateexnode = organ_finalstatepath.split('.')[0] + '.exnode'
        csvToExnode = CsvToExnode(organ_finalstatepath, organ_finalstateexnode, nodeOffset,
                                  self._organ_name + '_marker',
                                  coordinates_name=self._organ_name + '_marker_coordinates')
        csvToExnode.csv_to_exnode()
        self.read_landmarks(organ_finalstatepath)
        # self.get_element_xi(bodyNodeFile, bodyElemFile, organ_finalstateexnode)
        self.get_element_xi(bodyNodeFile, bodyElemFile, organ_finalstateexnode)

    def transform_organ(self, outputpath, *body_filepath):
        marker_list = self.get_marker_list()
        self.get_marker_coordinates(organ_referencepath, organ_filepath, marker_list=marker_list)
        # exelem = os.path.join(parentDirectory, 'rat1_fit1_element.exf')
        # exnode = os.path.join(parentDirectory, 'bodycoordAligned_node4.exf')
        # target_filename = os.path.join(working_directory, 'trans_{}0.csv'.format(organ_name))
        embedorgan.get_marker_coordinates(outputpath, *body_filepath, marker_list=marker_list)
        # findCoordInOtherTimes2(exnode, exelem, organ_name, parentDirectory, working_directory, target_filename, finalstatePath)
        target_filename = outputpath
        embedorgan.compute_transformation_matrix(working_directory, organ_name + '_reference.csv', target_filename)

    def get_marker_list(self):
        marker_list = []
        for landmark in self._organ_landmarks:
            marker_list.append(landmark['landmark'])

        return marker_list

# organ_name = 'brainstem'
# organ_filepath = r'C:\Users\egha355\Desktop\embed\human_brainstem\HumanBrainstemFit.exf'
# working_directory = r'C:\Users\egha355\Desktop\embed\human_brainstem\brainstem'
# parentDirectory = r'C:\Users\egha355\Desktop\embed\human_brainstem'
# bodyNodeFilepath = os.path.join(parentDirectory, 'human1_fit1_node1.exf')
# bodyElemFilepath = os.path.join(parentDirectory, 'human1_fit1_element.exf')
# #
# # body_nodefilepath = os.path.join(parentDirectory, 'rat1_fit1_node1.exf')
# finalstatefilename = organ_name + '_finalstate.csv'
# reference_filename = organ_name + '_reference.csv'
# # organ_finalstate = os.path.join(working_directory, organ_name + '_finalstate.exf')
# # organ_finalstateexnode = os.path.join(working_directory, organ_name + '_finalstate.exnode')
# organ_finalstatepath = os.path.join(working_directory, finalstatefilename)
# organ_referencepath = os.path.join(working_directory, reference_filename)
# dirname = os.path.dirname(organ_filepath)
# marker_filename = 'human1_fit1_node_static.exf'
# marker_filepath = os.path.join(parentDirectory, marker_filename)
#
# try:
#     organ_dirname = os.path.join(dirname, organ_name)
#     os.mkdir(organ_dirname)
# except OSError as error:
#     print(error)

organ_name = 'heart'
organ_filepath = r'C:\Users\egha355\Desktop\embed\embedding_static\mouse-with-organs\others\heart_transformed.exf'
working_directory = r'C:\Users\egha355\Desktop\embed\embedding_static\mouse-with-organs\others\heart'
parentDirectory = r'C:\Users\egha355\Desktop\embed\embedding_static\mouse-with-organs\others'
bodyNodeFilepath = os.path.join(parentDirectory, 'whole-body_node.exf')
bodyElemFilepath = os.path.join(parentDirectory, 'whole-body_element.exf')
#
# body_nodefilepath = os.path.join(parentDirectory, 'rat1_fit1_node1.exf')
finalstatefilename = organ_name + '_finalstate.csv'
reference_filename = organ_name + '_reference.csv'
# organ_finalstate = os.path.join(working_directory, organ_name + '_finalstate.exf')
# organ_finalstateexnode = os.path.join(working_directory, organ_name + '_finalstate.exnode')
organ_finalstatepath = os.path.join(working_directory, finalstatefilename)
organ_referencepath = os.path.join(working_directory, reference_filename)
dirname = os.path.dirname(organ_filepath)
marker_filename = 'human1_fit1_node_static.exf'
marker_filepath = os.path.join(parentDirectory, marker_filename)

try:
    organ_dirname = os.path.join(dirname, organ_name)
    os.mkdir(organ_dirname)
except OSError as error:
    print(error)

embedorgan = EmbedOrgan(organ_name, organ_filepath, organ_dirname, bodyNodeFilepath)
# Get the organ_finalstate.csv from fitted organ. This is only for fitted ones. For other cases the file should be given.
# embedorgan.get_marker_coordinates(organ_finalstatepath, organ_filepath)

embedorgan.embed_organ(organ_finalstatepath, organ_filepath, bodyNodeFilepath, bodyElemFilepath)
# embedorgan.get_marker_coordinates(finalstatePath, organ_filename)
# nodeOffset = embedorgan.get_maximum_node_identifier() + 100
# csvToExnode = CsvToExnode(finalstatePath, organ_finalstateexnode, nodeOffset, organ_name + '_marker', coordinates_name=organ_name+'_marker_coordinates')
# csvToExnode.csv_to_exnode()

# embedorgan.read_landmarks(organ_finalstatepath)
# embedorgan.get_element_xi(bodyNodeFile, bodyElemFile, organ_finalstateexnode)
# Copy element_xi to static file where we have marker and marker_location
# embedorgan.get_marker_coordinates(organ_referencepath, organ_filepath)
#
# marker_list = []
# for landmark in embedorgan._organ_landmarks:
#     marker_list.append(landmark['landmark'])

exelem = os.path.join(parentDirectory, 'human1_fit1_element.exf')
exnode = os.path.join(parentDirectory, 'align_node2.exf')
outputpath = os.path.join(working_directory, 'trans_{}0.csv'.format(organ_name))
embedorgan.transform_organ(outputpath, exnode, exelem, marker_filepath)
# embedorgan.get_marker_coordinates(target_filename, exnode, exelem, markerFilepath, marker_list=marker_list)
# findCoordInOtherTimes2(exnode, exelem, organ_name, parentDirectory, working_directory, target_filename, finalstatePath)
# embedorgan.compute_transformation_matrix(working_directory, organ_name + '_reference.csv', target_filename)

exnode = os.path.join(parentDirectory, 'scaffold_node3.exf')
outputpath = os.path.join(working_directory, 'trans_{}1.csv'.format(organ_name))
embedorgan.transform_organ(outputpath, exnode, exelem, marker_filepath)
# embedorgan.get_marker_coordinates(target_filename, exnode, exelem, markerFilepath, marker_list=marker_list)
# embedorgan.compute_transformation_matrix(working_directory, organ_name + '_reference.csv', target_filename)

exnode = os.path.join(parentDirectory, 'human1_fit1_node1.exf')
outputpath = os.path.join(working_directory, 'trans_{}2.csv'.format(organ_name))
embedorgan.transform_organ(outputpath, exnode, exelem, marker_filepath)
# embedorgan.get_marker_coordinates(target_filename, exnode, exelem, marker_filepath, marker_list=marker_list)
# embedorgan.compute_transformation_matrix(working_directory, organ_name + '_reference.csv', target_filename)

# exnode = os.path.join(parentDirectory, 'align_node2.exf')
# outputpath = os.path.join(working_directory, 'trans_{}3.csv'.format(organ_name))
# embedorgan.transform_organ(outputpath, exnode, exelem, marker_filepath)
# embedorgan.get_marker_coordinates(target_filename, exnode, exelem, marker_filepath, marker_list=marker_list)
# embedorgan.compute_transformation_matrix(working_directory, organ_name + '_reference.csv', target_filename)


text:
from pyspark.sql import SparkSession
from pyspark.sql.types import *
from pyspark.ml.feature import PCA
from pyspark.ml.clustering import KMeans
from pyspark.ml.linalg import Vectors
from pyspark.ml.feature import VectorAssembler
from pyspark.sql.functions import monotonically_increasing_id
from pyspark.sql import Row
from pyspark.ml.classification import MultilayerPerceptronClassifier
import numpy as np
import csv

# Create the spark object, so that the code can run on the spark.
spark = SparkSession \
    .builder \
    .appName("Python Spark KNN l3n200") \
    .getOrCreate() \

spark.conf.set("spark.sql.crossJoin.enabled", "true")

# define the path of data
# test_datafile = "hdfs://soit-hdp-pro-1.ucc.usyd.edu.au/share/MNIST/Test-label-28x28.csv"
# train_datafile= "hdfs://soit-hdp-pro-1.ucc.usyd.edu.au/share/MNIST/Train-label-28x28.csv"

# local path
test_datafile = "hdfs://soit-hdp-pro-1.ucc.usyd.edu.au/user/zzha7690/test_set.csv"
train_datafile = "hdfs://soit-hdp-pro-1.ucc.usyd.edu.au/user/zzha7690/train_set.csv"

# read the test_data_set and train_data_set
test_df = spark.read.csv(test_datafile, header=False, inferSchema="true")
train_df = spark.read.csv(train_datafile, header=False, inferSchema="true")

# transfer the test_df in to the dataframe with 2 column label and features so that we can do the further process
assembler_test = VectorAssembler(inputCols = test_df.columns[:1024],outputCol="fectures")
test_vectors_withlabel = assembler_test.transform(test_df).selectExpr("_c1024 as label_test","fectures")

# transfer the train_df in to the dataframe with 2 column label and features so that we can do the further process
assembler_train = VectorAssembler(inputCols = train_df.columns[:1024],outputCol="fectures")
train_vectors_withlabel = assembler_train.transform(train_df).selectExpr("_c1024 as label_train","fectures")

# fit the pca of train_vector first
# we set the k=50, so that we can keep 90% data of MINST
# After fit process, we can get the model of pca_200.
# Therefore, we can use the model to transform the test and train data.
pca = PCA(k=200, inputCol="fectures", outputCol="pca_vector")
model_200 = pca.fit(train_vectors_withlabel)
pca_train_result = model_200.transform(train_vectors_withlabel).selectExpr('label_train as label', 'pca_vector as feature')
pca_test_result = model_200.transform(test_vectors_withlabel).selectExpr('label_test as label', 'pca_vector as feature')

# define parameters
input_layer = 200 # number of features
output_layer = 10 # output 0~9
hidden_1 = 150
hidden_2 = 150
hidden_3 = 100
hidden_4 = 256
layers = [input_layer,hidden_1,hidden_1,output_layer]

MPC = MultilayerPerceptronClassifier(featuresCol='feature', labelCol='label', predictionCol='prediction', maxIter=100, layers=layers, blockSize=128, seed=123 )

model = MPC.fit(pca_train_result)

result = model.transform(pca_test_result).select("label", "prediction")
result_lp = result.selectExpr("label", "cast (prediction as int) prediction")
final_result = result_lp.rdd

# calculate the accuracy

neutral_zero_value = 0

def seqOp(a, b):
    if b[0] == b[1]:
        return a
    else:
        return a + 1

combOp = (lambda a, b: a + b)

diff = final_result.aggregate(neutral_zero_value, seqOp, combOp)
diff = float(diff)

accuracy = (26032 - diff) / 26032
print("The accuracy is", accuracy)

# Confusion matrix
confusion_matrix = np.zeros(shape=(10, 10))
confusion_matrix = confusion_matrix.astype(int)

def matrix_addvalue(matrix, value):
    matrix[int(value[1])][int(value[0])] = matrix[int(value[1])][int(value[0])] + 1
    return matrix

combine_matrix = (lambda m1, m2: m1 + m2)

result_confusion = final_result.aggregate(confusion_matrix, matrix_addvalue, combine_matrix)
with open('/home/zzha7690/ML_A2/Confusion_matrix.csv', 'w', newline='') as csvfile:
    spamwriter = csv.writer(csvfile, delimiter=',', quotechar='\n', quoting=csv.QUOTE_MINIMAL)
    spamwriter.writerows(result_confusion)

# Statistic table

def statistic(confusion_test):
    re_list = []
    label = np.arange(0, 10)
    for i in label:
        TP = confusion_test[i, i]
        FN = np.sum(confusion_test[i]) - TP
        FP = np.sum(confusion_test[:, i]) - TP
        TN = np.sum(confusion_test) - TP - FN - FP
        precision = (TP / (TP + FP))
        recall = TP / (TP + FN)
        F_measure = TP / (2 * TP + FP + FN)
        Support = (TP + FN)
        row = [int(label[i]), round(float(precision), 3), round(float(recall), 3), round(float(F_measure), 3),
               round(float(Support), 0)]
        re_list.append(row)
    return re_list

statistic_list = statistic(result_confusion)
statistic_rdd = spark.sparkContext.parallelize(statistic_list)
schema = StructType([
    StructField("Label", IntegerType(), True),
    StructField("Precision", FloatType(), True),
    StructField("recall", FloatType(), True),
    StructField("F1-score", FloatType(), True),
    StructField("Support", FloatType(), True)
])
statistic_df = spark.createDataFrame(statistic_rdd, schema)
statistic_df.show()

ext: py

reason: import numpy


text:
#!/usr/bin/python3

# generate .wav files containing per-class probabilities from .tf log files

# classify2.py <logdir> <model> <tffile> <context> <shift> <stride> [<labels>, <prevalences>]

# e.g.
# classify2.py `pwd`/trained-classifier 1k `pwd`/groundtruth-data/20161207T102314_ch1_p1.tf 204.8 0 1.6 mel-pulse,mel-sine,ambient 0.1,0.1,0.8

import sys
import os
import re
import numpy as np
from sys import argv
import ast
from scipy.io import wavfile

_,logdir,model,tffile,context_ms,shift_ms,stride_ms = argv[:7]
print('logdir: '+logdir)
print('model: '+model)
print('tffile: '+tffile)
print('context_ms: '+context_ms)
print('shift_ms: '+shift_ms)
print('stride_ms: '+stride_ms)
context_ms = float(context_ms)
shift_ms = float(shift_ms)
stride_ms = float(stride_ms)

with open(os.path.join(logdir,model,'vgg_labels.txt'), 'r') as fid:
  model_labels = fid.read().splitlines()

if len(argv)==9:
  labels,prevalences = argv[7:9]
  labels = np.array(labels.split(','))
  prevalences = np.array([float(x) for x in prevalences.split(',')])
  assert len(labels)==len(prevalences)==len(model_labels)
  assert set(labels)==set(model_labels)
  prevalences /= np.sum(prevalences)
  iimodel_labels = np.argsort(np.argsort(model_labels))
  ilabels = np.argsort(labels)
  labels = labels[ilabels][iimodel_labels]
  prevalences = prevalences[ilabels][iimodel_labels]
  assert np.all(labels==model_labels)
else:
  labels = model_labels
  prevalences = [1/len(labels) for _ in range(len(labels))]
print('labels: '+str(labels))
print('prevalences: '+str(prevalences))

npadding = round((context_ms/2+shift_ms)/stride_ms)
probability_matrix = np.zeros((npadding, len(labels)))
regex=re.compile('INFO:tensorflow:[0-9]+ms output_layer:0 (.+)')

with open(tffile) as fid:
  for line in fid:
    m = regex.match(line)
    if m==None:
      print(line)
      continue
    probabilities = ast.literal_eval(m.group(1).replace('nan','0'))
    probability_matrix = np.concatenate((probability_matrix,np.array(probabilities, ndmin=2)))

tic_rate = round(1000/stride_ms)
if tic_rate != 1000/stride_ms:
  print('WARNING: .wav files do not support fractional sampling rates!')

denominator = np.sum(probability_matrix * prevalences, axis=1)
for ch in range(len(labels)):
  adjusted_probability = probability_matrix[:,ch] * prevalences[ch]
  adjusted_probability[npadding:] /= denominator[npadding:]
  #inotnan = ~isnan(probability_matrix(:,ch));
  waveform = adjusted_probability*np.iinfo(np.int16).max  # colon was inotnan
  filename = tffile[:-3]+'-'+labels[ch]+'.wav'
  wavfile.write(filename, tic_rate, waveform.astype('int16'))

max_stars_repo_name: dstern/SongExplorer

max_issues_repo_name: dstern/SongExplorer

max_forks_repo_name: dstern/SongExplorer

reason: import numpy,from scipy


text:
import pytest
from numpy import allclose, arange, array, asarray, dot, cov, corrcoef, float64

from thunder.series.readers import fromlist, fromarray
from thunder.images.readers import fromlist as img_fromlist

pytestmark = pytest.mark.usefixtures("eng")

def test_map(eng):
    data = fromlist([array([1, 2]), array([3, 4])], engine=eng)
    assert allclose(data.map(lambda x: x + 1).toarray(), [[2, 3], [4, 5]])
    assert data.map(lambda x: 1.0*x, dtype=float64).dtype == float64
    assert data.map(lambda x: 1.0*x).dtype == float64

def test_map_singletons(eng):
    data = fromlist([array([4, 5, 6, 7]), array([8, 9, 10, 11])], engine=eng)
    mapped = data.map(lambda x: x.mean())
    assert mapped.shape == (2, 1)

def test_filter(eng):
    data = fromlist([array([1, 2]), array([3, 4])], engine=eng)
    assert allclose(data.filter(lambda x: x.sum() > 3).toarray(), [3, 4])

def test_flatten(eng):
    arr = arange(2*2*5).reshape(2, 2, 5)
    data = fromarray(arr, engine=eng)
    assert data.flatten().shape == (4, 5)
    assert allclose(data.flatten().toarray(), arr.reshape(2*2, 5))

def test_sample(eng):
    data = fromlist([array([1, 5]), array([1, 10]), array([1, 15])], engine=eng)
    assert allclose(data.sample(3).shape, (3, 2))
    assert allclose(data.filter(lambda x: x.max() > 10).sample(1).toarray(), [1, 15])

def test_between(eng):
    data = fromlist([array([4, 5, 6, 7]), array([8, 9, 10, 11])], engine=eng)
    val = data.between(0, 2)
    assert allclose(val.index, array([0, 1]))
    assert allclose(val.toarray(), array([[4, 5], [8, 9]]))

def test_first(eng):
    data = fromlist([array([4, 5, 6, 7]), array([8, 9, 10, 11])], engine=eng)
    assert allclose(data.first(), [4, 5, 6, 7])

def test_select(eng):
    index = ['label1', 'label2', 'label3', 'label4']
    data = fromlist([array([4, 5, 6, 7]), array([8, 9, 10, 11])], engine=eng, index=index)
    assert data.select('label1').shape == (2, 1)
    assert allclose(data.select('label1').toarray(), [4, 8])
    assert allclose(data.select(['label1']).toarray(), [4, 8])
    assert allclose(data.select(['label1', 'label2']).toarray(), array([[4, 5], [8, 9]]))
    assert data.select('label1').index == ['label1']
    assert data.select(['label1']).index == ['label1']

def test_standardize_axis1(eng):
    data = fromlist([array([1, 2, 3, 4, 5])], engine=eng)
    centered = data.center(1)
    standardized = data.standardize(1)
    zscored = data.zscore(1)
    assert allclose(centered.toarray(), array([-2, -1, 0, 1, 2]), atol=1e-3)
    assert allclose(standardized.toarray(),
                    array([0.70710,  1.41421,  2.12132,  2.82842,  3.53553]), atol=1e-3)
    assert allclose(zscored.toarray(),
                    array([-1.41421, -0.70710,  0,  0.70710,  1.41421]), atol=1e-3)

def test_standardize_axis0(eng):
    data = fromlist([array([1, 2]), array([3, 4])], engine=eng)
    centered = data.center(0)
    standardized = data.standardize(0)
    zscored = data.zscore(0)
    assert allclose(centered.toarray(), array([[-1, -1], [1, 1]]), atol=1e-3)
    assert allclose(standardized.toarray(), array([[1, 2], [3, 4]]), atol=1e-3)
    assert allclose(zscored.toarray(), array([[-1, -1], [1, 1]]), atol=1e-3)

def test_squelch(eng):
    data = fromlist([array([1, 2]), array([3, 4])], engine=eng)
    squelched = data.squelch(5)
    assert allclose(squelched.toarray(), [[0, 0], [0, 0]])
    squelched = data.squelch(3)
    assert allclose(squelched.toarray(), [[0, 0], [3, 4]])
    squelched = data.squelch(1)
    assert allclose(squelched.toarray(), [[1, 2], [3, 4]])

def test_correlate(eng):
    data = fromlist([array([1, 2, 3, 4, 5])], engine=eng)
    sig = [4, 5, 6, 7, 8]
    corr = data.correlate(sig).toarray()
    assert allclose(corr, 1)
    sigs = [[4, 5, 6, 7, 8], [8, 7, 6, 5, 4]]
    corrs = data.correlate(sigs).toarray()
    assert allclose(corrs, [1, -1])

def test_correlate_multiindex(eng):
    index = [[0, 1, 2, 3, 4], [5, 6, 7, 8, 9]]
    data = fromlist([array([1, 2, 3, 4, 5])], index=asarray(index).T, engine=eng)
    sig = [4, 5, 6, 7, 8]
    corr = data.correlate(sig).toarray()
    assert allclose(corr, 1)
    sigs = [[4, 5, 6, 7, 8], [8, 7, 6, 5, 4]]
    corrs = data.correlate(sigs).toarray()
    assert allclose(corrs, [1, -1])

def test_clip(eng):
    data = fromlist([array([1, 2, 3, 4, 5])], engine=eng)
    assert allclose(data.clip(2).toarray(), [2, 2, 3, 4, 5])
    assert allclose(data.clip(2, 3).toarray(), [2, 2, 3, 3, 3])

def test_mean(eng):
    data = fromlist([arange(8), arange(8)], engine=eng)
    val = data.mean().toarray()
    expected = data.toarray().mean(axis=0)
    assert allclose(val, expected)
    assert str(val.dtype) == 'float64'

def test_sum(eng):
    data = fromlist([arange(8), arange(8)], engine=eng)
    val = data.sum().toarray()
    expected = data.toarray().sum(axis=0)
    assert allclose(val, expected)
    assert str(val.dtype) == 'int64'

def test_var(eng):
    data = fromlist([arange(8), arange(8)], engine=eng)
    val = data.var().toarray()
    expected = data.toarray().var(axis=0)
    assert allclose(val, expected)
    assert str(val.dtype) == 'float64'

def test_std(eng):
    data = fromlist([arange(8), arange(8)], engine=eng)
    val = data.std().toarray()
    expected = data.toarray().std(axis=0)
    assert allclose(val, expected)
    assert str(val.dtype) == 'float64'

def test_max(eng):
    data = fromlist([arange(8), arange(8)], engine=eng)
    val = data.max().toarray()
    expected = data.toarray().max(axis=0)
    assert allclose(val, expected)

def test_min(eng):
    data = fromlist([arange(8), arange(8)], engine=eng)
    val = data.min().toarray()
    expected = data.toarray().min(axis=0)
    assert allclose(val, expected)

def test_labels(eng):
    x = [array([0, 1]), array([2, 3]), array([4, 5]), array([6, 7])]
    data = fromlist(x, labels=[0, 1, 2, 3], engine=eng)

    assert allclose(data.filter(lambda x: x[0]>2).labels, array([2, 3]))
    assert allclose(data[2:].labels, array([2, 3]))
    assert allclose(data[1].labels, array([1]))
    assert allclose(data[1, :].labels, array([1]))
    assert allclose(data[[0, 2]].labels, array([0, 2]))
    assert allclose(data.flatten().labels, array([0, 1, 2, 3]))

    x = [array([[0, 1],[2, 3]]), array([[4, 5], [6, 7]])]
    data = img_fromlist(x, engine=eng).toseries()
    data.labels = [[0, 1], [2, 3]]

    assert allclose(data.filter(lambda x: x[0]>1).labels, array([2, 3]))
    assert allclose(data[0].labels, array([[0, 1]]))
    assert allclose(data[:, 0].labels, array([[0], [2]]))
    assert allclose(data.flatten().labels, array([0, 1, 2, 3]))

def test_labels_setting(eng):
    x = [array([0, 1]), array([2, 3]), array([4, 5]), array([6, 7])]
    data = fromlist(x, engine=eng)

    with pytest.raises(ValueError):
        data.labels = [0, 1, 2]

def test_index_setting(eng):
    data = fromlist([array([1, 2, 3]), array([2, 2, 4]), array([4, 2, 1])], engine=eng)
    assert allclose(data.index, array([0, 1, 2]))
    data.index = [3, 2, 1]
    assert allclose(data.index, [3, 2, 1])
    with pytest.raises(ValueError):
        data.index = 5
    with pytest.raises(ValueError):
        data.index = [1, 2]

def test_select_by_index(eng):
    data = fromlist([arange(12)], index=[0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2], engine=eng)
    result = data.select_by_index(1)
    assert allclose(result.toarray(), array([4, 5, 6, 7]))
    assert allclose(result.index, array([1, 1, 1, 1]))
    result = data.select_by_index(1, squeeze=True)
    assert allclose(result.index, array([0, 1, 2, 3]))
    index = [
        [0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1],
        [0, 0, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1],
        [0, 1, 0, 1, 2, 3, 0, 1, 0, 1, 2, 3]
    ]
    data.index = array(index).T
    result, mask = data.select_by_index(0, level=2, return_mask=True)
    assert allclose(result.toarray(), array([0, 2, 6, 8]))
    assert allclose(result.index, array([[0, 0, 0], [0, 1, 0], [1, 0, 0], [1, 1, 0]]))
    assert allclose(mask, array([1, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0]))
    result = data.select_by_index(0, level=2, squeeze=True)
    assert allclose(result.toarray(), array([0, 2, 6, 8]))
    assert allclose(result.index, array([[0, 0], [0, 1], [1, 0], [1, 1]]))
    result = data.select_by_index([1, 0], level=[0, 1])
    assert allclose(result.toarray(), array([6, 7]))
    assert allclose(result.index, array([[1, 0, 0], [1, 0, 1]]))
    result = data.select_by_index(val=[0, [2,3]], level=[0, 2])
    assert allclose(result.toarray(), array([4, 5]))
    assert allclose(result.index, array([[0, 1, 2], [0, 1, 3]]))
    result = data.select_by_index(1, level=1, filter=True)
    assert allclose(result.toarray(), array([0, 1, 6, 7]))
    assert allclose(result.index, array([[0, 0, 0], [0, 0, 1], [1, 0, 0], [1, 0, 1]]))

def test_aggregate_by_index(eng):
    data = fromlist([arange(12)], index=[0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2], engine=eng)
    result = data.aggregate_by_index(sum)
    assert allclose(result.toarray(), array([6, 22, 38]))
    assert allclose(result.index, array([0, 1, 2]))
    index = [
        [0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1],
        [0, 0, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1],
        [0, 1, 0, 1, 2, 3, 0, 1, 0, 1, 2, 3]
    ]
    data.index = array(index).T
    result = data.aggregate_by_index(sum, level=[0, 1])
    assert allclose(result.toarray(), array([1, 14, 13, 38]))
    assert allclose(result.index, array([[0, 0], [0, 1], [1, 0], [1, 1]]))

def test_stat_by_index(eng):
    data = fromlist([arange(12)], index=[0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2], engine=eng)
    assert allclose(data.stat_by_index('sum').toarray(), array([6, 22, 38]))
    assert allclose(data.stat_by_index('mean').toarray(), array([1.5, 5.5, 9.5]))
    assert allclose(data.stat_by_index('min').toarray(), array([0, 4, 8]))
    assert allclose(data.stat_by_index('max').toarray(), array([3, 7, 11]))
    assert allclose(data.stat_by_index('count').toarray(), array([4, 4, 4]))
    assert allclose(data.stat_by_index('median').toarray(), array([1.5, 5.5, 9.5]))
    assert allclose(data.sum_by_index().toarray(), array([6, 22, 38]))
    assert allclose(data.mean_by_index().toarray(), array([1.5, 5.5, 9.5]))
    assert allclose(data.min_by_index().toarray(), array([0, 4, 8]))
    assert allclose(data.max_by_index().toarray(), array([3, 7, 11]))
    assert allclose(data.count_by_index().toarray(), array([4, 4, 4]))
    assert allclose(data.median_by_index().toarray(), array([1.5, 5.5, 9.5]))

def test_stat_by_index_multi(eng):
    index = [
        [0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1],
        [0, 0, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1],
        [0, 1, 0, 1, 2, 3, 0, 1, 0, 1, 2, 3]
    ]
    data = fromlist([arange(12)], index=array(index).T, engine=eng)
    result = data.stat_by_index('sum', level=[0, 1])
    assert allclose(result.toarray(), array([1, 14, 13, 38]))
    assert allclose(result.index, array([[0, 0], [0, 1], [1, 0], [1, 1]]))
    result = data.sum_by_index(level=[0, 1])
    assert allclose(result.toarray(), array([1, 14, 13, 38]))
    assert allclose(result.index, array([[0, 0], [0, 1], [1, 0], [1, 1]]))

def test_mean_by_panel(eng):
    data = fromlist([arange(8)], engine=eng)
    test1 = data.mean_by_panel(4)
    assert allclose(test1.index, array([0, 1, 2, 3]))
    assert allclose(test1.toarray(), [[2, 3, 4, 5]])
    test2 = data.mean_by_panel(2)
    assert allclose(test2.index, array([0, 1]))
    assert allclose(test2.toarray(), [[3, 4]])

def test_times_array(eng):
    mat1raw = asarray([[1, 2, 3], [4, 5, 6]])
    mat2 = asarray([[7, 8], [9, 10], [11, 12]])
    mat1 = fromlist(mat1raw, engine=eng)
    truth = dot(mat1raw, mat2)
    result = mat1.times(mat2)
    assert allclose(result.toarray(), truth)
    assert allclose(result.index, range(0, 2))

def test_times_array_alt(eng):
    mat1raw = asarray([[1, 2, 3], [4, 5, 6]])
    mat2 = asarray([[7, 8, 7, 8], [9, 10, 9, 10], [11, 12, 11, 12]])
    mat1 = fromlist(mat1raw, engine=eng)
    truth = dot(mat1raw, mat2)
    result = mat1.times(mat2)
    assert allclose(result.toarray(), truth)
    assert allclose(result.index, range(0, 4))

def test_times_vector(eng):
    mat1raw = asarray([[1, 2, 3], [4, 5, 6]])
    mat2 = [7, 8, 9]
    mat1 = fromlist(mat1raw, engine=eng)
    truth = dot(mat1raw, mat2)
    result = mat1.times(mat2)
    assert allclose(result.toarray(), truth)
    assert allclose(result.index, [0])

def test_times_scalar(eng):
    mat1raw = asarray([[1, 2, 3], [4, 5, 6]])
    mat2 = 5
    mat1 = fromlist(mat1raw, engine=eng)
    truth = mat1raw * mat2
    result = mat1.times(mat2)
    assert allclose(result.toarray(), truth)
    assert allclose(result.index, range(0, 3))

def test_gramian(eng):
    mat1raw = asarray([[1, 2, 3], [4, 5, 6]])
    mat1 = fromlist(mat1raw, engine=eng)
    result = mat1.gramian()
    truth = dot(mat1raw.T, mat1raw)
    assert allclose(result.toarray(), truth)

def test_cov(eng):
    mat1raw = asarray([[1, 2, 3], [4, 5, 6]])
    mat1 = fromlist(mat1raw, engine=eng)
    result = mat1.cov()
    truth = cov(mat1raw.T)
    assert allclose(result.toarray(), truth)

def test_fourier(eng):
    data = fromlist([array([1.0, 2.0, -4.0, 5.0, 8.0, 3.0, 4.1, 0.9, 2.3])], engine=eng)
    vals = data.fourier(freq=2)
    assert allclose(vals.select('coherence').toarray(), 0.578664)
    assert allclose(vals.select('phase').toarray(), 4.102501)

def test_convolve(eng):
    data = fromlist([array([1, 2, 3, 4, 5])], engine=eng)
    sig = array([1, 2, 3])
    betas = data.convolve(sig, mode='same')
    assert allclose(betas.toarray(), array([4, 10, 16, 22, 22]))

def test_crosscorr(eng):
    local = array([1.0, 2.0, -4.0, 5.0, 8.0, 3.0, 4.1, 0.9, 2.3])
    data = fromlist([local], engine=eng)
    sig = array([1.5, 2.1, -4.2, 5.6, 8.1, 3.9, 4.2, 0.3, 2.1])
    betas = data.crosscorr(signal=sig, lag=0)
    assert allclose(betas.toarray(), corrcoef(local, sig)[0, 1])
    betas = data.crosscorr(signal=sig, lag=2)
    truth = array([-0.18511, 0.03817, 0.99221, 0.06567, -0.25750])
    assert allclose(betas.toarray(), truth, atol=1E-5)

def test_detrend(eng):
    data = fromlist([array([1, 2, 3, 4, 5])], engine=eng)
    out = data.detrend('linear')
    assert(allclose(out.toarray(), array([1, 1, 1, 1, 1])))

def test_normalize_percentile(eng):
    data = fromlist([array([1, 2, 3, 4, 5])], engine=eng)
    out = data.normalize('percentile', perc=20)
    vals = out.toarray()
    assert str(vals.dtype) == 'float64'
    assert allclose(vals, array([-0.42105,  0.10526,  0.63157,  1.15789,  1.68421]), atol=1e-3)

def test_normalize_window(eng):
    y = array([1, 2, 3, 4, 5])
    data = fromlist([y], engine=eng)
    vals = data.normalize('window', window=2).toarray()
    b = array([1, 1, 2, 3, 4])
    result_true = (y - b) / (b + 0.1)
    assert allclose(vals, result_true, atol=1e-3)
    vals = data.normalize('window', window=5).toarray()
    b = array([1, 1, 2, 3, 4])
    result_true = (y - b) / (b + 0.1)
    assert allclose(vals, result_true, atol=1e-3)

def test_normalize_window_exact(eng):
    y = array([1, 2, 3, 4, 5])
    data = fromlist([y], engine=eng)
    vals = data.normalize('window-exact', window=2).toarray()
    b = array([1.2,  1.4,  2.4,  3.4,  4.2])
    result_true = (y - b) / (b + 0.1)
    assert allclose(vals, result_true, atol=1e-3)
    vals = data.normalize('window-exact', window=6).toarray()
    b = array([1.6,  1.8,  1.8,  1.8,  2.6])
    result_true = (y - b) / (b + 0.1)
    assert allclose(vals, result_true, atol=1e-3)

def test_normalize_mean(eng):
    data = fromlist([array([1, 2, 3, 4, 5])], engine=eng)
    vals = data.normalize('mean').toarray()
    assert allclose(vals, array([-0.64516,  -0.32258,  0.0,  0.32258,  0.64516]), atol=1e-3)

def test_mean_by_window(eng):
    data = fromlist([array([0, 1, 2, 3, 4, 5, 6])], engine=eng)
    test1 = data.mean_by_window(indices=[3, 5], window=2).toarray()
    assert allclose(test1, [3, 4])
    test2 = data.mean_by_window(indices=[3, 5], window=3).toarray()
    assert allclose(test2, [3, 4, 5])
    test3 = data.mean_by_window(indices=[3, 5], window=4).toarray()
    assert allclose(test3, [2, 3, 4, 5])
    test4 = data.mean_by_window(indices=[3], window=4).toarray()
    assert allclose(test4, [1, 2, 3, 4])

def test_reshape(eng):
    original =  fromarray(arange(72).reshape(6, 6, 2), engine=eng)
    arr = original.toarray()

    assert allclose(arr.reshape(12, 3, 2), original.reshape(12, 3, 2).toarray())
    assert allclose(arr.reshape(36, 2), original.reshape(36, 2).toarray())
    assert allclose(arr.reshape(4, 3, 3, 2), original.reshape(4, 3, 3, 2).toarray())

    # must conserve number of elements
    with pytest.raises(ValueError):
        original.reshape(6, 3, 2)

    # cannot change length of series
    with pytest.raises(ValueError):
        original.reshape(6, 3, 4)

def test_downsample(eng):
    data = fromlist([arange(8)], engine=eng)
    vals = data.downsample(2).toarray()
    assert allclose(vals, [0.5, 2.5, 4.5, 6.5])
    vals = data.downsample(4).toarray()
    assert allclose(vals, [1.5, 5.5])

def test_downsample_uneven(eng):
    data = fromlist([arange(9)], engine=eng)
    vals = data.downsample(2).toarray()
    assert allclose(vals, [0.5, 2.5, 4.5, 6.5])

max_issues_repo_name: gopikasula/thunder

max_forks_repo_name: gopikasula/thunder

reason: from numpy


text:
import random

import numpy as np

def metropolis(func, rv, n, downsample=1):
    """
    Metropolis algorithm

    Parameters
    ----------
    func : callable
        (un)normalized distribution to be sampled from
    rv : RandomVariable
        proposal distribution which is symmetric at the origin
    n : int
        number of samples to draw
    downsample : int
        downsampling factor

    Returns
    -------
    sample : (n, ndim) ndarray
        generated sample
    """
    x = np.zeros((1, rv.ndim))
    sample = []
    for i in range(n * downsample):
        x_new = x + rv.draw()
        accept_proba = func(x_new) / func(x)
        if random.random() < accept_proba:
            x = x_new
        if i % downsample == 0:
            sample.append(x[0])
    sample = np.asarray(sample)
    assert sample.shape == (n, rv.ndim), sample.shape
    return sample

max_stars_repo_name: alexandru-dinu/PRML

max_issues_repo_name: alexandru-dinu/PRML

max_forks_repo_name: alexandru-dinu/PRML


text:
# -*- coding: utf-8 -*-
"""
Created on Thu Feb 18 16:26:24 2021

@author: hayua
"""
import matplotlib.pyplot as plt

from pyswarm import pso
"Imports"
import CoolProp.CoolProp as CP
import numpy as np

from misc.utils import NoStdStreams
from components.heat_exchangers import hx

def model(x):

    P_ev=x[0]
    P_cd=x[1]
    p = np.zeros(8, dtype=float)
    v = np.zeros(8, dtype=float)
    T = np.zeros(8, dtype=float)
    x = np.zeros(8, dtype=float)
    h = np.zeros(8, dtype=float)
    s = np.zeros(8, dtype=float)

    # Cycle parameters:
    fluid ='propane'
    # P_ev = 6000000  #pa
    # P_cd = 200000  #pa
    p_low=P_cd
    p_high=P_ev

    DELTAT_sh = 5
    DELTAT_sc = 5
    epsilon_s = 0.6

    # Heat sink parameters:
    T_hsin_in =273.15-160           # Natural gas inlet temperature /K
    DELTAT_cd=10                    # pinch point temperature of cooling water /K
    # Q_dot_cd=10000                  # heat capacity flowrates of cooling water kW/K

    # Heat source parameters:
    T_hsor_in=273.15 + 200        # cooling water inlet temperature /K
    DELTAT_pinch=10                    # pinch point temperature of cooling water /K

    # get fluid properties:
    #T_crit = CP.PropsSI("Tcrit",fluid)
    #p_crit = CP.PropsSI("Pcrit",fluid)
    T_low = CP.PropsSI('T','Q', 0.5, 'T', P_ev, fluid)
    T_high = CP.PropsSI('T','Q', 0.5, 'T', P_cd, fluid)
    #h_crit = CP.PropsSI('H','T', T_crit, 'P', p_crit, fluid)
    #s_crit = CP.PropsSI('S','T', T_crit, 'P', p_crit, fluid)

    #Evaporator outlet:
    p[0] = p_low
    s[0] = CP.PropsSI('S','Q', 0, 'P', p_low, fluid)
    T[0] = CP.PropsSI('T','Q', 0, 'P', p_low, fluid)
    h[0] = CP.PropsSI('H','Q', 0, 'P', p_low, fluid)

    # pump inlet:
    p_su_cp = p_low
    T_su_cp = CP.PropsSI('T','Q', 0, 'P', p_low, fluid)-DELTAT_sc
    h_su_cp = CP.PropsSI('H','T', T_su_cp, 'P', p_su_cp, fluid)
    s_su_cp = CP.PropsSI('S','T', T_su_cp, 'P', p_su_cp, fluid)
    s[1] = s_su_cp
    T[1] = T_su_cp
    h[1] = h_su_cp
    p[1] = p_su_cp

    #pump outlet:
    p_ex_cp = p_high
    h_ex_cp_s = CP.PropsSI('H','S', s_su_cp, 'P', p_ex_cp, fluid)
    h_ex_cp = h_su_cp - (h_su_cp - h_ex_cp_s)/epsilon_s
    T_ex_cp = CP.PropsSI('T','H', h_ex_cp, 'P', p_ex_cp, fluid)
    s_ex_cp = CP.PropsSI('S','H', h_ex_cp, 'P', p_ex_cp, fluid)

    p[2] = p_high
    s[2] = s_ex_cp
    T[2] = T_ex_cp
    h[2] = h_ex_cp

    #From statepoint 2 to state ponit 3, we assume seawater is used to heat to state point 3 at 0 Celsius:
    p[3] = p_high
    T[3] = 273.15
    s[3] = CP.PropsSI('S','T', T[3] , 'P', p_high, fluid)
    h[3] = CP.PropsSI('H','T', T[3] , 'P', p_high, fluid)

    p[4] = p_high
    T[4]=190+273.15
    s[4] = CP.PropsSI('S','T', T[4] , 'P', p_high, fluid)
    h[4] = CP.PropsSI('H','T', T[4] , 'P', p_high, fluid)

    p[5] = p_low
    h_ex_turb_s = CP.PropsSI('H','S', s[4], 'P', p_low, fluid)
    h[5] = h[4] - (h[4] - h_ex_turb_s)*epsilon_s
    T[5] = CP.PropsSI('T','H', h[5], 'P', p_low, fluid)
    s[5] = CP.PropsSI('S','H', h[5], 'P', p_low, fluid)

    p[6] = p_low
    s[6] = CP.PropsSI('S','Q', 1, 'P', p_low, fluid)
    T[6] = CP.PropsSI('T','Q', 1, 'P', p_low, fluid)
    h[6] = CP.PropsSI('H','Q', 1, 'P', p_low, fluid)

    p[7] = p[0]
    s[7] = s[0]
    T[7] = T[0]
    h[7] = h[0]

    # heat sink:
    cp_w = 100000  # Heat capactity flow rate of both hot water and LNG 100kW/C
    M_dot1 = cp_w*(T[6]-DELTAT_pinch-T_hsin_in)/ (h[6] - h[1]) # mass flowrate calculated by LNG end

    M_dot=M_dot1
    W_tur=(h[4]-h[5])*M_dot
    W_pump=(h[2]-h[1])*M_dot
    W_net=W_tur-W_pump

    heatld_eva=M_dot*(h[4]-h[3])
    T_hsor_out=T_hsor_in-heatld_eva/cp_w
    heatld_cond=M_dot*(h[5]-h[1])
    T_hsin_out=T_hsin_in+heatld_cond/cp_w
    heatld_cond_sup=M_dot*(h[5]-h[6])

    # print("The temperature of each state")
    # print(T)
    # print("The pressure of each state (Pa) ")
    # print(p)
    # print("The enthalpy of each state (J/kg)")
    # print(h)
    # print("The entropy of each state (J/kg/K)")
    # print(s)
    # print ("The mass flow rate of working fluid is %.2f kg/s" %(M_dot))
    # print("The condensation temperature is %.2f C" %(T[6]-273.15))
    # print("The heat load of evaproator is %.2f kW and the outlet temperature of heat source is %.2f C" %(heatld_eva/1000,T_hsor_out-273.15))
    # print("The heat load of condenser is %.2f kW and the outlet temperature of heat sink is %.2f C" %(heatld_cond/1000,T_hsin_out-273.15))
    # print("The turbine power output is %.2f kW" %(W_tur/1000))
    # print("The pump power input is %.2f kW" %(W_pump/1000))
    # print("The net power output is %.2f kW" %(W_net/1000))
    # Heat source:
    # print(heatld_cond_sup)

    # Temperature profile in the condenser:
    # print("For the condenser:")
    s_cd,T_cd,T_hf,pinch_cd = hx(fluid,100,M_dot,h[5],h[1],p_low,p_low,T_hsin_in,cp_w)
    # plt.plot(T_cd)
    # plt.plot(T_hf)
    # plt.show()
    # # # Temperature profile in the evaporator
    # print("For the evaporator:")
    s_ev,T_ev,T_cf,pinch_ev = hx(fluid,100,M_dot,h[3],h[4],p_high,p_high,T_hsor_in,cp_w)
    # plt.plot(T_ev)
    # plt.plot(T_cf)
    # plt.show()

    while pinch_ev<10: # if the pinch temperature is less than 5, we just decrease the target temperature of working fluid

        T[4]-=2   #step size is chosen as 2 celsius
        s[4] = CP.PropsSI('S','T', T[4] , 'P', p_high, fluid)
        h[4] = CP.PropsSI('H','T', T[4] , 'P', p_high, fluid)
        h_ex_turb_s = CP.PropsSI('H','S', s[4], 'P', p_low, fluid)
        h[5] = h[4] - (h[4] - h_ex_turb_s)*epsilon_s
        T[5] = CP.PropsSI('T','H', h[5], 'P', p_low, fluid)
        s[5] = CP.PropsSI('S','H', h[5], 'P', p_low, fluid)

        M_dot=M_dot1
        W_tur=(h[4]-h[5])*M_dot
        W_pump=(h[2]-h[1])*M_dot
        W_net=W_tur-W_pump

        heatld_eva=M_dot*(h[4]-h[3])
        T_hsor_out=T_hsor_in-heatld_eva/cp_w
        heatld_cond=M_dot*(h[5]-h[1])
        T_hsin_out=T_hsin_in+heatld_cond/cp_w
        heatld_cond_sup=M_dot*(h[5]-h[6])
        s_cd,T_cd,T_hf,pinch_cd = hx(fluid,100,M_dot,h[5],h[1],p_low,p_low,T_hsin_in,cp_w)
        s_ev,T_ev,T_cf,pinch_ev = hx(fluid,100,M_dot,h[3],h[4],p_high,p_high,T_hsor_in,cp_w)

    print("The temperature of each state")
    print(T)
    print("The pressure of each state (Pa) ")
    print(p)
    print("The enthalpy of each state (J/kg)")
    print(h)
    print("The entropy of each state (J/kg/K)")
    print(s)
    print ("The mass flow rate of working fluid is %.2f kg/s" %(M_dot))
    print("The condensation temperature is %.2f C" %(T[6]-273.15))
    print("The heat load of evaproator is %.2f kW and the outlet temperature of heat source is %.2f C" %(heatld_eva/1000,T_hsor_out-273.15))
    print("The heat load of condenser is %.2f kW and the outlet temperature of heat sink is %.2f C" %(heatld_cond/1000,T_hsin_out-273.15))
    print("The turbine power output is %.2f kW" %(W_tur/1000))
    print("The pump power input is %.2f kW" %(W_pump/1000))
    print("The net power output is %.2f kW" %(W_net/1000))

    print("For the condenser:")
    s_cd,T_cd,T_hf,pinch_cd = hx(fluid,100,M_dot,h[5],h[1],p_low,p_low,T_hsin_in,cp_w)
    plt.plot(T_cd)
    plt.plot(T_hf)
    plt.show()
    # # Temperature profile in the evaporator
    print("For the evaporator:")
    s_ev,T_ev,T_cf,pinch_ev = hx(fluid,100,M_dot,h[3],h[4],p_high,p_high,T_hsor_in,cp_w)
    plt.plot(T_ev)
    plt.plot(T_cf)
    plt.show()

    D=CP.PropsSI('D','T|gas',T[5],'P',p[5],fluid)
    V=M_dot/D
    v=V/W_net

    from CoolProp.Plots import PropertyPlot
    from CoolProp.Plots.SimpleCycles import StateContainer
    import pickle
    import os

    cache_plot = False
    if cache_plot:
        filename = fluid + '.p'
        if os.path.isfile(filename):   # Load previously saved plot
            pp = pickle.load(open( filename, "rb" ) )
        else:
            states = StateContainer()
            states_hf = StateContainer()
            states_cf = StateContainer()
            pp = PropertyPlot('HEOS::'+fluid, 'HEOS::R134a','TS', tp_limits='DEF')
            pp.set_axis_limits([-0.5, 3, -200, 530])
            with NoStdStreams():
                pp.calc_isolines()
            with open(filename, 'wb') as f:
                pickle.dump(pp, f)
    else:
        states = StateContainer()
        states_hf = StateContainer()
        states_cf = StateContainer()
        pp = PropertyPlot('HEOS::'+fluid, 'TS')
        pp.set_axis_limits([0, 3.5, 100, 530])
        with NoStdStreams():
            pp.calc_isolines()

    for i in range(4):
        states[i,'T'] = T[i]
        states[i,"S"] = s[i]

    for i,Tx in enumerate(T_ev):
        states.append({'T':Tx,'S':s_ev[i]})

    for i in range(4,len(T)):
        states.append({'T':T[i],'S':s[i]})
    states.append({'T':T[1],'S':s[1]})    # for some reasons, the second point needs to be repeated to close the cycle

    for i,Tx in enumerate(T_hf):
        states_hf.append({'T':Tx,'S':s_cd[i]})

    for i,Tx in enumerate(T_cf):
        states_cf.append({'T':Tx,'S':s_ev[i]})

    with NoStdStreams():
        pp.draw_process(states,line_opts={'color':'green'})
        pp.draw_process(states_hf,line_opts={'color':'blue', 'linestyle':'dashed'})
        pp.draw_process(states_cf,line_opts={'color':'red', 'linestyle':'dashed'})

    pp.show()

    return -W_net,V,v

lb = [8129000, 650000] #lower bound with CO2 as the working fluid
ub = [15000000,3000000]  #upper bound with CO2 as the working fluid

# lb = [4675000, 100000] #lower bound with propane as the working fluid
# ub = [10000000,1500000]  #upper bound with propane as the working fluid

# xopt, fopt = pso(model, lb, ub,swarmsize=50, maxiter=100)
# xopt=[11569607.3265, 1161490.312] # optimum value with evaporation pressure upper bound 150 bar with CO2 as the working fluid

xopt=[5863765.83, 100000]  #optimum value with evaporation pressure upper bound 150 bar with propane as the working fluid

netpower,Volume_total, Volume_specific=model(xopt)

# xopt=[10000000, 998484] currently obtained best results.
#%%

# from CoolProp.Plots import PropertyPlot
# from CoolProp.Plots.SimpleCycles import StateContainer
# import pickle
# import os

# cache_plot = False
# if cache_plot:
#     filename = fluid + '.p'
#     if os.path.isfile(filename):   # Load previously saved plot
#         pp = pickle.load(open( filename, "rb" ) )
#     else:
#         states = StateContainer()
#         states_hf = StateContainer()
#         states_cf = StateContainer()
#         pp = PropertyPlot('HEOS::'+fluid, 'TS')
#         with NoStdStreams():
#             pp.calc_isolines()
#         with open(filename, 'wb') as f:
#             pickle.dump(pp, f)
# else:
#     states = StateContainer()
#     states_hf = StateContainer()
#     states_cf = StateContainer()
#     pp = PropertyPlot('HEOS::'+fluid, 'TS')
#     with NoStdStreams():
#         pp.calc_isolines()

# for i in range(3):
#     states[i,'T'] = T[i]
#     states[i,"S"] = s[i]

# for i,Tx in enumerate(T_cd):
#     states.append({'T':Tx,'S':s_cd[i]})

# for i in range(4,len(T)):
#     states.append({'T':T[i],'S':s[i]})
# states.append({'T':T[1],'S':s[1]})    # for some reasons, the second point needs to be repeated to close the cycle

# for i,Tx in enumerate(T_hf):
#     states_hf.append({'T':Tx,'S':s_cd[i]})

# for i,Tx in enumerate(T_cf):
#     states_cf.append({'T':Tx,'S':s_ev[i]})

# with NoStdStreams():
#     pp.draw_process(states,line_opts={'color':'green'})
#     pp.draw_process(states_hf,line_opts={'color':'red', 'linestyle':'dashed'})
#     pp.draw_process(states_cf,line_opts={'color':'blue', 'linestyle':'dashed'})

# pp.show()

max_stars_repo_path: Transcritical power cycle.py

max_stars_repo_name: squoilin/LNG-Regasification

max_issues_repo_path: Transcritical power cycle.py

max_issues_repo_name: squoilin/LNG-Regasification

max_forks_repo_path: Transcritical power cycle.py

max_forks_repo_name: squoilin/LNG-Regasification


text:
#!/usr/bin/env python
import numpy as np
import pandas as pd
from sklearn.utils import check_random_state
from .base import TreatmentDescriptionDataFrame, \
    DelayedDataFrame
from ..utils._validation import _check_dframe, _check_feature_target, \
    _check_train_test_size, _check_column_dtypes, \
    _check_boolean, \
    _check_int_keyword, _check_float_keyword, _check_str_keyword, \
    _check_min_feature_significance, _check_hc_threshold, \
    _check_ordinal_mapping, _check_cat_fill_value

class TreatmentDesignMixin(object):
    """
    Mixin class defining core functionality for treatment design-related
    classes.
    """

    def __init__(self, features="all", target=-1,
                 min_feature_significance=None, cv=3, cv_type='loo',
                 cv_split_function=None, train_size=0.2,
                 test_size=None, random_state=None,
                 unique_level_min_percent=0.02,
                 novel_level_strategy="nan", rare_level_pooling=False,
                 make_nan_indicators=True,
                 high_cardinality_strategy="impact",
                 downstream_context=None, remove_duplicates=False,
                 feature_scaling=False,
                 rare_level_significance=None,
                 feature_engineering=True,
                 convert_dtypes=True,
                 find_ordinals=True,
                 ordinals=None,
                 high_cardinality_threshold=None,
                 n_jobs=1,
                 ordinal_mapping=None,
                 categorical_fill_value="NaN",
                 ordinal_fill_value=-1):

        # TODO: ensure all params, attrs, etc. are documented
        self.features = features
        self.target = target
        self.significance = min_feature_significance
        self.cv = cv
        self.cv_type = cv_type
        self.split_function = cv_split_function
        self.train_size = train_size
        self.test_size = test_size
        self.random_state = random_state
        self.unique_percent = unique_level_min_percent
        self.novel_strategy = novel_level_strategy
        self.make_nans = make_nan_indicators
        self.rare_pooling = rare_level_pooling
        self.high_cardinality_strategy = high_cardinality_strategy
        self.high_cardinality_threshold = high_cardinality_threshold
        self.context = downstream_context
        self.deduplicate = remove_duplicates
        self.scaling = feature_scaling
        self.rare_significance = rare_level_significance
        self.feature_engineering = feature_engineering
        self.convert_dtypes = convert_dtypes
        self.find_ordinals = find_ordinals
        self.ordinals=ordinals
        self.n_jobs = n_jobs
        self.mapping = ordinal_mapping
        self.categorical_fill_value = categorical_fill_value
        self.ordinal_fill_value = ordinal_fill_value

    def fit(self, dataframe):
        """
        Create the treatment plan.
        """
        # input validation
        self._validate_params(dataframe)

        # with input validated, create the TreatmentDescriptionDF object
        self.TreatmentDescriptionDF_ = TreatmentDescriptionDataFrame()
        self.DelayedDF_ = DelayedDataFrame()

        # copy the dataframe for modification
        self.df_ = dataframe.copy()

        # get the indices for nan columns
        self._get_nan_indices()

        # fill in missing values
        self._fill_na()

        # self._get_categorical_value_counts()
        # make indicators for additional modeling and significance testing
        self._make_indicators()

        # make indicators for whether a value replacement occurred
        self._make_replacement_indicators(dataframe)

        return self

    def transform(self, dataframe):
        """
        Apply the treatment to the dataframe.
        """
        pass

    def fit_transform(self, dataframe):
        """
        Create and apply the treatment plan.
        """
        pass

    def _validate_params(self, dataframe):
        """
        Validates the parameters passed in during initialization.
        """
        _check_dframe(dataframe=dataframe)

        _check_min_feature_significance(self.significance,
                                        ["float between 0.0 and 1.0",
                                         "1/n_features",
                                         None])

        # TODO: add enhancement to allow value to be set based upon the dataset instead of fixed int values
        _check_hc_threshold(self.high_cardinality_threshold,
                            ["int >= 0", None])

        _check_ordinal_mapping(self.mapping,
                               [dict, None],
                               dataframe.columns.tolist())

        _check_cat_fill_value(self.categorical_fill_value,
                              "categorical_fill_value")

        _check_float_keyword(self.unique_percent,
                             "unique_level_min_percent",
                             [0, 1])

        _check_int_keyword(self.n_jobs,
                           "n_jobs",
                           [-1, np.inf])

        _check_int_keyword(self.ordinal_fill_value,
                           "ordinal_fill_value",
                           None)

        _check_str_keyword(self.context,
                           "downstream_context",
                           ["pipeline", None])

        _check_str_keyword(self.high_cardinality_strategy,
                           "high_cardinality_strategy",
                           ["impact", "indicators"])

        _check_str_keyword(self.novel_strategy,
                           "novel_level_strategy",
                           ["known", "zero", "nan", "rare", "pooled"])

        _check_boolean(self.rare_pooling, "rare_level_pooling")
        _check_boolean(self.make_nans, "make_nan_indicators")
        _check_boolean(self.find_ordinals, "find_ordinals")
        _check_boolean(self.deduplicate, "remove_duplicates")
        _check_boolean(self.scaling, "feature_scaling")
        _check_boolean(self.feature_engineering, "feature_engineering")

        self.remove_features_ = False if self.significance is None else True

        self.features_, self.target_ = _check_feature_target(
            df_columns=dataframe.columns,
            features=self.features,
            target=self.target)

        self.train_size_ = _check_train_test_size(self.train_size,
                                                  self.test_size)

        self.numeric_, self.categorical_, self.ordinal_, self.mapping_, \
        self.numeric_has_nan_, self.categorical_has_nan_, self.ordinal_has_nan_ = \
            _check_column_dtypes(dataframe=dataframe,
                                 features=self.features_,
                                 target=self.target_,
                                 convert_dtypes=self.convert_dtypes,
                                 ordinals=self.ordinals,
                                 find_ordinals=self.find_ordinals,
                                 categorical_fill_value=self.categorical_fill_value)

        self.random_state_ = check_random_state(self.random_state)

        # TODO: finish validation routines
        # TODO: check cv, cv_type, cv_split_function, rare_level_significance, convert_dtypes, ordinals

    def _get_nan_indices(self):
        # get the indices of missing values
        self.nan_indices_ = dict()

        self.nan_indices_["numeric_has_nan_"] = dict()
        self.nan_indices_["categorical_has_nan_"] = dict()
        self.nan_indices_["ordinal_has_nan_"] = dict()

        for column in self.numeric_has_nan_:
            self.nan_indices_["numeric_has_nan_"][column] = \
                np.where(self.df_.loc[:, column].isnull())[0]

        for column in self.categorical_has_nan_:
            self.nan_indices_["categorical_has_nan_"][column] = \
                np.where(self.df_.loc[:, column].isnull())[0]

        for column in self.ordinal_has_nan_:
            self.nan_indices_["ordinal_has_nan_"][column] = \
                np.where(self.df_.loc[:, column].isnull())[0]

    def _fill_na(self):
        # fill in missing values
        if len(self.numeric_has_nan_) > 0:
            self.nan_numeric_ = self.numeric_has_nan_ + "_clean"

            self.df_ = self.df_.reindex(columns=[self.df_.columns.tolist() +
                                                 self.nan_numeric_.tolist()],
                                        fill_value=0)

            for c1, c2 in zip(self.nan_numeric_, self.numeric_has_nan_):
                self.df_.loc[:, c1] = self.df_.loc[:, c2].copy()
                self.df_.loc[:, c1] = self.df_.loc[:, c1].fillna(
                    self.df_.loc[:, c1].mean())
        else:
            self.nan_numeric_ = pd.Index(list())

        # no need to drop these as we're turning them into indicators anyway
        # which will drop them
        if len(self.categorical_has_nan_) > 0:

            self.df_.loc[:, self.categorical_has_nan_] = \
                self.df_.loc[:, self.categorical_has_nan_] \
                    .fillna(value=self.categorical_fill_value)
        else:
            self.nan_categorical_ = None

        # TODO: apply mapping for ordinals if present
        # TODO: we should probably look earlier for mix/max ordinal values and set the fill value based upon that to prevent collisions
        # TODO: this will have to change if we start allowing ordinal detection for non-int dtypes
        if len(self.ordinal_has_nan_) > 0:
            self.nan_ordinal_ = self.ordinal_has_nan_ + "_clean"
            self.df_ = self.df_.reindex(columns=[self.df_.columns.tolist() +
                                                 self.nan_ordinal_.tolist()],
                                        fill_value=0)

            for c1, c2 in zip(self.nan_ordinal_, self.ordinal_has_nan_):
                self.df_.loc[:, c1] = self.df_.loc[:, c2].copy()
                self.df_.loc[:, c1] = self.df_.loc[:, c1].fillna(
                    self.df_.loc[:, c1].mean())
        else:
            self.nan_ordinal_ = None

    def _make_indicators(self):
        # make indicators for replaced values for all dtypes
        if self.categorical_ is None:
            self.nan_categorical_ = None
        else:
            self.nan_categorical_ = self.categorical_has_nan_ + "_" + self.categorical_fill_value
            self.df_ = pd.get_dummies(self.df_,
                                      columns=self.categorical_.tolist())

    def _make_replacement_indicators(self, dataframe):
        if self.nan_numeric_ is not None:
            self.numeric_replacement_indicators_ = self.numeric_has_nan_ + \
                                                   "_replaced"

            # add the rep indicator columns to the dataframe
            self.df_ = \
                self.df_.reindex(columns=[self.df_.columns.tolist() +
                                          self.numeric_replacement_indicators_.tolist()],
                                 fill_value=0)

            # add the indicator if that row value was replaced
            for c1, c2 in zip(self.numeric_replacement_indicators_,
                              self.numeric_has_nan_):
                self.df_.loc[:, c1] = dataframe.loc[:, c2]\
                    .isnull()\
                    .map({True: 1,
                          False: 0})

            # drop the original columns with NaN
            self.df_.drop(labels=self.numeric_has_nan_.tolist(),
                          inplace=True,
                          axis=1)

        if self.nan_categorical_ is not None:
            self.categorical_replacement_indicators_ = self.categorical_has_nan_ + \
                                                       "_replaced"

            # add the rep indicator columns to the dataframe
            self.df_ = \
                self.df_.reindex(columns=[self.df_.columns.tolist() +
                                          self.categorical_replacement_indicators_.tolist()],
                                 fill_value=0)

            # add the indicator if that row value was replaced
            for c1, c2 in zip(self.categorical_replacement_indicators_,
                              self.categorical_has_nan_):
                self.df_.loc[:, c1] = dataframe.loc[:, c2].isnull().map(
                    {True: 1,
                     False: 0})

        if self.nan_ordinal_ is not None:
            self.ordinal_replacement_indicators_ = self.ordinal_has_nan_ + \
                                                   "_replaced"

            # add the rep indicator columns to the dataframe
            self.df_ = \
                self.df_.reindex(columns=[self.df_.columns.tolist() +
                                          self.ordinal_replacement_indicators_.tolist()],
                                 fill_value=0)

            # add the indicator if that row value was replaced
            for c1, c2 in zip(self.ordinal_replacement_indicators_,
                              self.ordinal_has_nan_):
                self.df_.loc[:, c1] = dataframe.loc[:, c2].isnull().map(
                    {True: 1,
                     False: 0})

            # drop the original columns with NaN
            self.df_.drop(labels=self.ordinal_has_nan_.tolist(),
                          inplace=True,
                          axis=1)

max_stars_repo_name: jmwoloso/avalearn

max_issues_repo_name: jmwoloso/avalearn

max_forks_repo_name: jmwoloso/avalearn


text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
#   http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing,
# software distributed under the License is distributed on an
# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied.  See the License for the
# specific language governing permissions and limitations
# under the License.

# Authors of this page : Beatriz Kanzki & Victor Dupuy

#!/usr/bin/Python-2.7.11

from django.shortcuts import render
from django.http import HttpResponse
from django.http import JsonResponse
from django.db import connection # Used to connect with the database

#--- Pandas ---#
import pandas as pandas
import numpy as numpy

#--- Regex ---#
import re as regex

#--- CSV reader ---#
import csv

#--- Tools ---#
from Tools import connect
from Tools import getPhenotypes
from Tools import UserLogs

#--- JSON ---#
import json

import urllib2

def adamGenomeViewer(request, chromosome, position, rsid, userWidth, userHeight):
    print request, chromosome, position, rsid

    params = "rsid=" + rsid + "&chromosome=" + chromosome + "&position=" + position;

    data = urllib2.urlopen("http://ec2-52-35-68-107.us-west-2.compute.amazonaws.com:8000/matching?" + params).read()

    rsidArray = []
    rsidArray.append(json.loads(data))#[o['rsid'] for o in data]

    chrBoundaries = getChromosomeBoundaries()

    rsidArray = json.dumps(rsidArray, ensure_ascii=False, encoding="utf-8").replace("\\", "")

    response = json.dumps({
            'data' : {
                'jsonChrBoundaries' : chrBoundaries,
                'jsonValidRsids' : rsidArray
            }
        },
        sort_keys=True,
        indent=4,
        separators=(',', ': ')
    )

    return HttpResponse(response)

def sqlGenomeViewer(data):

    rsidArray = data

    chrBoundaries = getChromosomeBoundaries()
    validRsids = fetchValidRsids(rsidArray)

    print "validRsids"
    print validRsids

    response = json.dumps({
            'data' : {
                'jsonChrBoundaries' : chrBoundaries,
                'jsonValidRsids' : validRsids
            }
        },
        sort_keys=True,
        indent=4,
        separators=(',', ': ')
    )

    return HttpResponse(response)

def getChromosomeBoundaries():
    #We query the database for min and max positions for each chromosome
    sqlQuery = "SELECT chromosome, MIN(position) as min, MAX(position) as max FROM marqueurs GROUP BY chromosome;"
    chrBoundaries = connect.fetchData(sqlQuery)
    return buildJsonData(chrBoundaries)

def fetchValidRsids(rows):
    #This code was used back when we only searched for matching rsids (no position or chromosome).
    #print "rsids"
    #print rsids
    #print ','.join(map(str, rsids))
    #We query all rows where we have a data match
    #sqlQuery = 'SELECT * FROM marqueurs WHERE nom in ("' + '","'.join(map(str, rsids)) + '");'
    #validRsids = connect.fetchData(sqlQuery)

    validRsids = []

    for row in rows:
        sqlQuery = 'SELECT * FROM marqueurs WHERE nom="' + row['rsid'] + '" AND position=' + row['position'] + ' AND chromosome=' + row['chromosome'] + ';'
        match = connect.fetchData(sqlQuery)

        # We do [1:-1] because it extracts the data from stringified-array. The data being a string, it removes the first and last characters : [ and ].
        j = buildJsonData(match)[1:-1]
        # If no result was returned, the data was "[]" which both chars were stripped. So.. empty string.
        if j != "":
            validRsids.append(json.loads(j))

    print "validRsids"
    print validRsids

    return json.dumps(validRsids)

def buildJsonData(data):
    jsonData = data.to_json(orient='records')#json table
    return jsonData

def uploadFile(request):
    return HttpResponse()

def extractHeader(request):
    f = request.FILES['file']
    reader = csv.reader(f)
    headers = reader.next()
    response = json.dumps({'headers': headers})
    f.close()
    return HttpResponse(response)

def fileGenomeViewer(request):
    f = request.FILES['file']
    reader = csv.DictReader(f)
    data = json.dumps([ row for row in reader ])
    f.close()

    output = []

    for row in json.loads(data):
        output.append({
            'rsid': row[request.POST['rsid_header']],
            'position': row[request.POST['position_header']],
            'chromosome': row[request.POST['chromosome_header']]
        })

    return sqlGenomeViewer(output)


text:
#! usr/bin/env python
# coding:utf-8
#=====================================================
# Copyright (C) 2020 * Ltd. All rights reserved.
#
# Author      : Chen_Sheng19
# Editor      : VIM
# Create time : 2020-05-07
# File name   : train.py
# Description : for training model with VOC2007 dataset
#
#=====================================================

import os
import time
import shutil
import numpy as np
import tensorflow as tf
from tqdm import tqdm
import utils
import dataset
from yolo3 import yolo3,decode,loss_layer
import getConfig

cfg = {}
cfg = getConfig.getConfig()

trainset = dataset.Dataset('train')
logdir = "./model_dir/log"
steps_per_epoch = len(trainset)
global_steps = tf.Variable(1,trainable = False,dtype=tf.int64)
warmup_steps = cfg['train_warmup_epochs']
total_steps = cfg['train_epochs'] * steps_per_epoch

input_tensor = tf.keras.layers.Input([416,416,3])
conv_tensors = yolo3(input_tensor)

output_tensor = []
for i,conv_tensor in enumerate(conv_tensors):
    pred_tensor = decode(conv_tensor,i)
    output_tensor.append(conv_tensor)
    output_tensor.append(pred_tensor)

model = tf.keras.Model(input_tensor,output_tensor)

utils.load_partial_weights(model,"./yolov3.weights")

for layer in model.layers[0:58]:
    layer.trainable = False

for layer in model.layers[59:66]:
    layer.trainable = False

for layer in model.layers[67:74]:
    layer.trainable = False
model.summary()

opt = tf.keras.optimizers.Adam()

if os.path.exists(logdir):
    shutil.rmtree(logdir)
writer = tf.summary.create_file_writer(logdir)

def train_step(image_data,target):
    with tf.GradientTape() as tape:
        pred_result = model(image_data)
        giou_loss = conf_loss = prob_loss = 0

        for i in range(3):
            conv,pred = pred_result[i*2],pred_result[i*2+1]
            loss_items = loss_layer(conv,pred,target[i],target[i+3],i)
            giou_loss += loss_items[0]
            conf_loss += loss_items[1]
            prob_loss += loss_items[2]

        total_loss = giou_loss + conf_loss + prob_loss

        gradients = tape.gradient(total_loss,model.trainable_variables)
        opt.apply_gradients(zip(gradients,model.trainable_variables))

        tf.print("=> STEP %4d  Learning-rate: %.6f giou-loss: %4.2f conf-loss: %4.2f prob-loss: %4.2f total-loss: %4.2f"%(global_steps,
            opt.lr.numpy(),giou_loss,conf_loss,prob_loss,total_loss))

        global_steps.assign_add(1)
        if global_steps < warmup_steps:
            lr = global_steps / warmup_steps * cfg['train_learning_rate_init']
        else:
            lr = cfg['train_learning_rate_end'] + 0.5 *(cfg['train_learning_rate_init'] - cfg['train_learning_rate_end']) * (
                    1 + tf.cos((global_steps - warmup_steps) / (total_steps - warmup_steps) * np.pi))
        opt.lr.assign(lr.numpy())

        with writer.as_default():
            tf.summary.scalar("lr",opt.lr,step = global_steps)
            tf.summary.scalar("loss/total_loss",total_loss,step = global_steps)
            tf.summary.scalar("loss/giou_loss",giou_loss,step = global_steps)
            tf.summary.scalar("loss/conf_loss",conf_loss,step = global_steps)
            tf.summary.scalar("loss/prob_loss",prob_loss,step = global_steps)
        writer.flush()

for epoch in range(cfg['train_epochs']):
    for image_data,*target in trainset:
        train_step(image_data,target)
    model.save_weights("./model_dir/yolo3")