text:
'''
This program is designed to gather to generate the required Pathset.xml file.
Primarily by collecting the .mts files that are output by sheba when measuring SWS, for each phase a seperate <data> tag is required
Stating the paths to the requisit SAC files (which are copied to the data directory in the current model path [for now at least])
#
'''
######## Imports ##############
from xml.etree import ElementTree
# ElementTree is a standard (as of python 2.5) library which we can use to parse XML.
# N.B ET is NOT secure against "malicous XML" however as we are only intersted in very simple XML this shouldnt be an issue
from xml.dom import minidom
import pandas as pd
import os
import glob
import numpy as np
from sphe_trig import vincenty_dist
from calc_aoi import slw2aoi,get_rayparam
from sactools import get_sac,get_mts
from corrections import add_sside_correction, add_station_correction,bin2domain
###############################
__author__ = "Joseph Asplet"
__license__ = "MIT"
__version__ = "0.1"
__email__ = "joseph.asplet@bristol.ac.uk"
__status__ = "Development"

class PathSetter:
    """A class to hold the metadata for the run (rdir, station? [for now], outdir etc. ) and fucntions
    to parse/generate the XML needed for the PathSet file
    """
    def __init__(self,phasefile, domains=None,model=None,odir=None,config_uid='Test Run'):
        '''
        Initialise Pathsetter with some essential data, metadata.

        Args:
            phasefile (str) - file that contains splitting results for all phases
            phases (list) - list of phases that are to be used
            domains (str) - file that contains information on the size/ shape of the domains
            model (str) - model XML file.If the model file does not exists Pathsetter can create one
            odir (str) - the output directory that XML files are written to
            config_uid (str) - a unique identifier written to MTSConfig.xml, default is 'Test Run'

        Returns:
            Set (obj) - a configured Pathsetter Object.

        Examples:
            >>> Set = Pathsetter('E_pacific_05_binned.pairs',['SKS','SKKS'],
                                'ScS_w_bouncepoints_binned_05.sdb','ScS','TestModel.xml')
        '''

        print('Reading df')
        self.date_time_convert = {'TIME': lambda x: str(x),'DATE': lambda x : str(x)}
        self.df = pd.read_csv(phasefile,converters=self.date_time_convert,delim_whitespace=True)
        print('Read Trigonal Domians file')
        if domains:
            self.doms = pd.read_csv(domains,delim_whitespace=True).sort_values(by='BIN')
        else:
            print('Not using trigonal bins')
        print('Set Outdir')
        if odir == None:
            self.opath = os.getcwd() # Gets current working directory and stores is as a Path
            self.odir = self.opath.split('/')[-1]
        elif len(odir.split('/')) > 1:
            self.opath = odir
            self.odir = odir.split('/')[-1]
        else:
            self.opath = '/Users/ja17375/SWSTomo/BlueCrystal/{}'.format(odir)
            self.odir = odir
        # Set XML namespace
        self.xmlns = {'mtsML':'http://www1.gly.bris.ac.uk/cetsei/xml/MatisseML/'}
        if model == None:
            print('Model will need to be generated before Pathsetting can be done')
        else:
            print("Reading Model file from SWSTomo/Models")
            self.modelxml = '/Users/ja17375/SWSTomo/Models/{}'.format(model)
            self.modelroot = ElementTree.parse('{}'.format(self.modelxml)).getroot()
            self.model = self.modelroot.find('mtsML:model',self.xmlns)
            model_name = self.model[0].text
            print('Using model named ... {}'.format(model_name))
        # Set config uID
        self.config_uid = config_uid

    def domain2operator(self,domain,phase,evdp,gcarc,azi):
        '''Function to read model.xml and extract the required domain and generate the operator XML needed for the Pathset file

        Args:
            domain (str) - the domain uid for the domain to make the operator for
            phase (str) - the phase passing through the domain
            evdp (str) - source event depth for phase
            gcarc (float) - event-station distance [deg]
            azi (float) - event station azimuth
        Returns:
            operator (obj) - An ElementTree Strucutre conainting the operator XML
        '''
        domains = self.model.findall('mtsML:domain',self.xmlns)
        uid = None
        for dom in domains:
            # Loop through domains in the model file
#             uid_tmp = dom.find('mtsML:domain_uid',self.xmlns).text
            uid_tmp = dom[0].text
            if domain == uid_tmp:
                uid = uid_tmp
                break

        if uid is None:
            raise ValueError('Domain {} Not Found'.format(domain))

        if uid.split('_')[0] == 'Lower':
            # Domains starting with Lower are at CMB. So depth == 2890 km
            depth = 2890. # Approx depth of CMB [km]

            slw = get_rayparam(evdp,gcarc,phase)
            aoi = slw2aoi(depth,slw) # Calculate ray param and then incidence angle
            dom_h = 250
            if phase == 'ScS':
                # double distance if ScS to account for the up- and -downgoing legs
                dist = (2*dom_h) / np.cos(np.deg2rad(aoi))
            else:
                dist = dom_h / np.cos(np.deg2rad(aoi))
        elif (uid.split('_')[0] in ['Upper','RSide','Station']):
            aoi = 0
            dist = 100 # 100km so we can use conversion dt/100 = strength
        elif uid.split('_')[0] == 'SSide':
            # If domain ID is SSide then we know this domain is for a source side correction.
            # azi = 0 # For the splitting corrections we fix azi to 0
            aoi = 0
            dist = 100
        else:
            raise ValueError(f'Domain "{uid}"name does not match expected values')

        operator = ElementTree.Element('operator')
        dom_uid = ElementTree.SubElement(operator,'domain_uid')
        dom_uid.text = domain
        azimuth = ElementTree.SubElement(operator,'azi')
        azimuth.text = str(azi)
        # azimuth.text = '0'
        inclination = ElementTree.SubElement(operator,'inc')
        inclination.text = str(90 - aoi) # change from aoi to inclination
        l = ElementTree.SubElement(operator,'dist')
        l.text = str(np.around(dist,decimals=3))
        return operator

    def parsedoms(self):
        '''
        Function to parse the input model XML into lists of domains for use in other functions

        Args:
            None - model XML is read in __init__ and is an attribute of the Setter Class (self)
        Returns:
            rdoms (list) - a list of all upper mantle reciever side domains in the model
            ldoms (list) - a list of all D`` domains in the model
            sdoms (list) - a list of all upper mantle source side domains in the model
        Examples:
            >>>(rdoms, ldoms, sdoms) = Set.parsedoms()
        '''
        dml = '{{{}}}domain'.format(self.xmlns['mtsML'])
        # need the triple curly braces to get a set of braces surrounding the mtsML
        rdoms = []
        ldoms = []
        sdoms = []
        self.udom_c = {}
        self.ldom_c = {}
        self.sdom_c = {}
        for dom in self.model.iter(dml):
            d_id = dom[0].text
            if d_id.split('_')[0] == 'Station':
                rdoms.append(d_id)
                self.udom_c.update({d_id : 0})
            elif d_id.split('_')[0] == 'Lower':
                ldoms.append(d_id)
                self.ldom_c.update({d_id : 0})
            elif d_id.split('_')[0] == 'SSide':
                # If the uid is SSide we know it is a source side correction domain. We just need to
                # match the correct domain to the correct phase
                sdoms.append(d_id)
                self.sdom_c.update({d_id : 0 })

        return (rdoms,ldoms,sdoms)

    def is_phase_good(self,row,split_c=0.7, null_c=-0.9):
        '''
        Function to test if a phase is not split or null based on its Q value (Wuestefeld et al., 2010)

        Args:
            row (df) - a row from a pandas DataFrame containing data for multiple phases
            ph (str) - the phase to do the Q test for

        Returns:
            res (bool) - The result of the Q test. True if phase passed, False if it fails.
        '''
        if (row['Q'] < 0.7) and (row['Q'] > -0.9):
            res = False
        else:
            res = True
        return res

    def loop_thru_domains(self,doms,raylat,raylon,uid,ph,crit,evdp,gcarc,azi):
        '''
        This function loops through a set of input domains and identifies domains whose midpoints are within a set distance from a raypath (the point the raypath is at at a given depth in the mantle). For these domains it then creates the required operator XML sequence

        Args:
            doms (list) - list of domains to loop over
            raylat (float) - latitude of the raypath
            raylon (float) - longitude of the raypath
            uid (str) - domain uid, for creating the XML (e.g. 'Upper')
            ph (str) - the phase ID ('SKS','SKKS','ScS')
            crit (float) - the distance ciriterea we want to find domains for
            evdp (str) - source event depth for phase
            gcarc (float) - event-station distance [deg]
            azi (float) - event station azimuth

        Returns:
            operators (list) - a list containing operator objects from domain2operator for the domains that were within the test distance from the raypath.
'''
        operators = [ ]
        for d,row in doms.iterrows():
            domlat = row.MID_LAT
            domlon = row.MID_LON
            dom_id = "{}_{:.0f}".format(uid,row.BIN)
            dist, azi = vincenty_dist(raylat,raylon,domlat,domlon)
            if dist <= crit:
                op = self.domain2operator(dom_id,ph,evdp,gcarc,azi)
                operators.append(op)

        return operators

    def gen_PathSet_XML(self,stations, lower_dom_no=None, phases=['ScS','SKS','SKKS'],fname=None, f=1):
        '''
        Function to iterate over the shear-wave phase data provided to Setter and to generate the correct XML
        to from the Pathset XML file required for Matisse. This includes adding source/ reciever side corrections
        where required (and where the neccesary corrections have been provided as needed by functions called
                        by gen_PathSet_XML)
        This funciton is written assuming there is a singular D'' domain to invert for

        Args:
            stations [list] - list of stations to include data for
            lower_dom_no [int or list] - lower domain of interest
            phases [list] - list of phase codes to iterate over. Default is ['ScS','SKS','SKKS']
            fname [str] - File name to write the Pathset XML to. If no name is provided then this defaults to 'Pathset'
        Returns:
            Output XML is written to fname
        '''
        # Some quick i/o management
        if fname == None:
            self.pathset_xml = 'Pathset'
        else:
            self.pathset_xml = fname

        #Start of main function
        self.pathset_root = ElementTree.Element('MatisseML')
        tree = ElementTree.ElementTree(self.pathset_root)
        self.pathset_root.set("xmlns",self.xmlns['mtsML'])
        pathset = ElementTree.SubElement(self.pathset_root,'pathset')
        psuid = 'Paths for run in dir {} .'.format(self.odir)
        pathset_uid = ElementTree.SubElement(pathset,'pathset_uid')
        pathset_uid.text = psuid
        p = 0
        # Before Pathsetting, parse the upper/lower domains from the model file
        rdoms,ldoms,sdoms = self.parsedoms()
        if stations == 'all':
            df = self.df
        elif stations == 'random_sample':
            df = self.df.sample(frac=f, replace=True)
        else:
            df = self.df[self.df.STAT.isin(stations)]

        for i, row in df.iterrows():
            # All XML generation must sit within this loop (function calls) so that we make sure that Az, EVDP etc. are right for the current phase
            attribs = {'evdp':row.EVDP,'azi':row.AZI,'gcarc':row.DIST,'evla':row.EVLA,
                            'evlo':row.EVLO,'stla':row.STLA,'stlo':row.STLO,
                            'date':row.DATE,'time':row.TIME,'stat':row.STAT,
                            'phase' :row.PHASE}
            filepath = self.set_filepaths(attribs)
            try:
                fileID = glob.glob(filepath)[0].strip('.mts').split('/')[-1]
                # Strip out .mts and split by '/', select end to get filestem
            except IndexError:
                continue
            statdom = 'Station_{}'.format(attribs['stat'])
            stat_op = self.domain2operator(statdom,attribs['phase'],attribs['evdp'],
                                               attribs['gcarc'],0)
            if attribs['phase'] == 'ScS': # Only ScS needs source side domains
                dom_id = 'SSide_{}_{}_{}'.format(attribs['stat'],
                                                 attribs['date'],attribs['time'])
                sside_op = self.domain2operator(dom_id,'ScS',attribs['evdp'],
                                                attribs['gcarc'],0)
            _, lmm_azi = vincenty_dist(row.LOWMM_LAT, row.LOWMM_LON, row.STLA, row.STLO)
            ldom_id = 'Lower_{}'.format(lower_dom_no)
            lowmm_op = self.domain2operator(ldom_id,attribs['phase'],attribs['evdp'],
                                   attribs['gcarc'],lmm_azi)

            get_sac(fileID,attribs['stat'],attribs['phase'])
            # Now make XML for this Path
            path = ElementTree.SubElement(pathset,'path')
            pathname = 'Path {} {}'.format(p,attribs['phase'])
            path_uid = ElementTree.SubElement(path,'path_uid')
            path_uid.text = pathname
            # Add Data (from .mts)
            data = get_mts(fileID,attribs['stat'],attribs['phase'])
            path.append(data)
            stat_uid = ElementTree.SubElement(path,'station_uid')
            stat_uid.text = attribs['stat']
            evt_uid = ElementTree.SubElement(path,'event_uid')
            evt_uid.text = '{}_{}'.format(attribs['date'],attribs['time'])
            if attribs['phase'] == 'ScS':
                # Add SSide Operator to path
                path.append(sside_op)

            # Add D`` Operator
            path.append(lowmm_op)
            path.append(stat_op)
            p+=1
        #Now write out the Pathset XML
        self._write_pretty_xml(self.pathset_root,file='{}/{}.xml'.format(self.opath,self.pathset_xml))

    def gen_synth_PathSet(self,lower_dom_no=None, fname=None, syntest='RealBAZ/Noise01'):
        '''
        Function to iterate over the shear-wave phase data provided to Setter and to generate the correct XML
        to from the Pathset XML file required for Matisse. This includes adding source/ reciever side corrections
        where required (and where the neccesary corrections have been provided as needed by functions called
                        by gen_PathSet_XML)
        This funciton is written assuming there is a singular D'' domain to invert for

        Args:
            stations [list] - list of stations to include data for
            lower_dom_no [int or list] - lower domain of interest
            phases [list] - list of phase codes to iterate over. Default is ['ScS','SKS','SKKS']
            fname [str] - File name to write the Pathset XML to. If no name is provided then this defaults to 'Pathset'
        Returns:
            Output XML is written to fname
        '''
        # Some quick i/o management
        if fname == None:
            self.pathset_xml = 'Pathset'
        else:
            self.pathset_xml = fname

        #Start of main function
        self.pathset_root = ElementTree.Element('MatisseML')
        tree = ElementTree.ElementTree(self.pathset_root)
        self.pathset_root.set("xmlns",self.xmlns['mtsML'])
        pathset = ElementTree.SubElement(self.pathset_root,'pathset')
        psuid = 'Paths for run in dir {} .'.format(self.odir)
        pathset_uid = ElementTree.SubElement(pathset,'pathset_uid')
        pathset_uid.text = psuid
                # Before Pathsetting, parse the upper/lower domains from the model file
        _ = self.parsedoms()
        for i, row in self.df.iterrows():
            # All XML generation must sit within this loop (function calls) so that we make sure that Az, EVDP etc. are right for the current phase
            attribs = {'evdp':row.EVDP,'azi':row.AZI,'gcarc':row.DIST,'evla':row.EVLA,
                            'evlo':row.EVLO,'stla':row.STLA,'stlo':row.STLO,
                            'date':row.DATE,'time':row.TIME,'stat':row.STAT,
                            'phase' :row.PHASE, 'station':row.STAT}
            filepath = f'/Users/ja17375/DiscrePy/Sheba/Runs/SYNTH/InversionSYNTH/{syntest}'
            try:
                test = '{}/{}*.mts'.format(filepath, attribs['station'])
                filestem = glob.glob(test)[0].strip('.mts').split('/')[-1]
                fileID = f'{filepath}/{filestem}'
                # Strip out .mts and split by '/', select end to get filestem
            except IndexError:
                continue
            ldom_id = 'Lower_{}'.format(lower_dom_no)
            lowmm_op = self.domain2operator(ldom_id,attribs['phase'],attribs['evdp'],
                                   attribs['gcarc'],attribs['azi'])

            get_sac(fileID,attribs['stat'],attribs['phase'], syn=True)
            # Now make XML for this Path
            path = ElementTree.SubElement(pathset,'path')
            pathname = 'Path Syn {}'.format(attribs['station'])
            path_uid = ElementTree.SubElement(path,'path_uid')
            path_uid.text = pathname
            # Add Data (from .mts)
            data = get_mts(fileID,attribs['station'],attribs['phase'], syn=True)
            path.append(data)
            stat_uid = ElementTree.SubElement(path,'station_uid')
            stat_uid.text = attribs['stat']
            evt_uid = ElementTree.SubElement(path,'event_uid')
            evt_uid.text = '{}_{}'.format(attribs['date'],attribs['time'])
            # Add D`` Operator
            path.append(lowmm_op)

        #Now write out the Pathset XML
        self._write_pretty_xml(self.pathset_root,file='{}/{}.xml'.format(self.opath,self.pathset_xml))

    def set_filepaths(self,attribs):
        '''Set filepaths to datasets for different phases
            I need this to keep up with the sprawling sheba Runs directories for different things
        '''
        if attribs['phase'] == 'ScS':
            f = '/Users/ja17375/DiscrePy/Sheba/Runs/ScS/{}/{}/{}_{}_{}??_{}.mts'.format(
                attribs['stat'],attribs['phase'],attribs['stat'],attribs['date'],attribs['time'],attribs['phase'])
        else:
            f = '/Users/ja17375/DiscrePy/Sheba/Runs/E_pacific/{}/{}/{}_{}_{}??_{}.mts'.format(attribs['stat'],
                attribs['phase'],attribs['stat'],attribs['date'],attribs['time'],attribs['phase'])
        return f

    def gen_Model_XML(self, stations, mod_name=None, Low_Domains=None, ):
        '''
        Generates the model XML file required by matisse, with lowermost mantle domains and seperate upper mantle domains for every station requested. Source Side corrections for ScS are auto-added

        Args:
            stations (lsit) - list of all Stations to create domains for.
            mod_name (str) - name of the model. The model XML will be written to a file mod_name.xml
            Low_Domains (array-like) - [optional] array of D`` domains to use in model. If none we assume this is a single domain run
        Returns:
            model XML written to the file mod_name.xml
        '''
        if mod_name is None:
            mod_name = input('No model name provided. Enter one now :')
        modpath = '/Users/ja17375/SWSTomo/Models'
        # Write the initial xml
        root = ElementTree.Element('MatisseML')
        tree = ElementTree.ElementTree(root)
        root.set("xmlns",self.xmlns['mtsML'])
        root.append(ElementTree.Comment('Generated by gen_Model from pathset.py'))
        m = ElementTree.SubElement(root,'model')
        m_uid = ElementTree.SubElement(m,'model_uid')
        m_uid.text = 'Trigonal Domains'

        if Low_Domains is None:
            print('No Lowermost Mantle Domains')

        if stations == 'all':
            stations =  self.df.STAT.values

        for stat in stations:
            dom = add_station_correction(stat)
            m.append(dom)

        for i,row in self.df.iterrows():
            if (row.STAT in stations) & (row.PHASE == 'ScS'):
                date, time, stat = row.DATE, row.TIME, row.STAT
                print(f'{date}, {time}, {stat}')
                sdom = add_sside_correction(date, time, stat)
                m.append(sdom)

        if Low_Domains:
            for ldom in Low_Domains:
                # Loop over requested lower domains
                print(ldom)
                dom = bin2domain(f'Lower_{ldom}')
                m.append(dom)

        print('Model Generated')
        self.model = m
        self._write_pretty_xml(root,file='{}/{}.xml'.format(modpath,mod_name))

    def count_paths_in_model(self,Low_Domains=None,Rside_Domains=None):
        '''
        This function takes a test set of D'' and Reciever Side UM domains (default is all available domains) and counts the number of ScS, SKS and SKKS phases that are within the criteria distance from each of them. This is so that we can get an idea of how many paths are sampling each domain before we generate all the XML. This function repeats some things done in gen_Pathset and gen_Model, but here we do not have to worry about any XML. All we want to output is an array of domains and the number of paths for each phase that pass through the domain

        Args:
            Low_Domains (array-like) - [optional] array of D`` domains to use in model. If none, all domains that have some phases in the input data to Setter are used (IF you have used geogeom to bin the data already)
            Rside_Domains (array-like) - [optional] array of all reciever-side domains to use in model. If none, reciever side for all phases (from the input data to the Setter class) assigned to a D`` domain will be added.

        Returns:
            counts (array-like) - a 2d numpy array (shape(1280,7)) that contains the counts for the number of phases that pass through each D`` and reciever side domain requested.
        '''
        date_time_convert = {'TIME': lambda x: str(x),'DATE': lambda x : str(x)}
        df = pd.read_csv('E_pacific_SNR10_goodQ.sdb',converters=date_time_convert,delim_whitespace=True)
        crit = 5.0 # [deg] the criteria distance
        if Low_Domains is None:
            ldoms = self.doms.BIN.values
        else:
            ldoms = np.asarray(Low_Domains)

        if Rside_Domains is None:
            rdoms = self.doms.BIN.values
        else:
            rdoms = np.asarray(Rside_Domains)

        doms = self.doms[self.doms.BIN.isin(ldoms) | self.doms.BIN.isin(rdoms)]
        print('Counting paths in domains {}'.format(doms.BIN.values))

        counts = np.zeros([1280,4])

        for i,dom in doms.iterrows():
            for i, row in df.iterrows():
                dlat, dlon = dom.MID_LAT, dom.MID_LON
                dom_id = int(dom.BIN)
                dist_lowmm = vincenty_dist(row.LOWMM_LAT, row.LOWMM_LON, dlat, dlon)[0]
#                 dist_skks = vincenty_dist(row.SKKS_LM_LAT, row.SKKS_LM_LON, dlat, dlon)[0]
                dist_rside = vincenty_dist(row.STLA, row.STLO, dlat, dlon)[0]
                if (dist_lowmm <= crit) & (dom.BIN in ldoms):
                    if row.PHASE == 'ScS':
                        counts[dom_id-1,0] = dom_id
                        counts[dom_id-1,1] += 1
                    elif (row.PHASE == 'SKS') | (row.PHASE == 'SKKS'):
                        counts[dom_id-1,0] = dom_id
                        counts[dom_id-1,2] += 1
                    else:
                        raise ValueError('Phase code {} incorrect'.format(row.phase))

        mask = np.all(np.equal(counts, 0), axis=1)
        c_out = counts[~mask]
        print('Done')
        return c_out

    def gen_MTS_Info(self):
        '''
        Function to create the simple XML file "MTS_Info"
        '''
        root = ElementTree.Element('MatisseML')
        tree = ElementTree.ElementTree(root)
        root.set("xmlns",self.xmlns['mtsML'])
        root.append(ElementTree.Comment('Generated by gen_MTS_Info from pathset.py'))
        samplerinfo = ElementTree.SubElement(root,'SamplerInfo')
        c_uid = ElementTree.SubElement(samplerinfo,'config_uid')
        c_uid.text = self.config_uid
        # Now write the XML to a file
        self._write_pretty_xml(root,file='{}/MTS_Info.xml'.format(self.opath))

    def gen_MTS_Config(self,options=None):
        '''
        Function to create the MTS config file. N.B this function will just generate a default config file. It is easier to edit/copy an existing file as only a few changes are needed to tweak a run.

        Args:
            options (dict) - options to be set in the Config file. N.B all options must be specified

        '''
        model = self.modelxml.split('/')[0]

        root = ElementTree.Element('MatisseML')
        tree = ElementTree.ElementTree(root)
        root.set("xmlns",self.xmlns['mtsML'])
        root.append(ElementTree.Comment('Generated by gen_MTS_Config from pathset.py'))
        config = ElementTree.SubElement(root,'config')
        c_uid = ElementTree.SubElement(config,'config_uid')
        c_uid.text = self.config_uid
        config.append(ElementTree.Comment(' input files'))
        mf = ElementTree.SubElement(config,'model_file') # Creat model file tag
        mf.text = model
        ps = ElementTree.SubElement(config,'pathset_file')
        ps.text = '{}.xml'.format(self.pathset_xml)
        config.append(ElementTree.Comment('operational options')) # Break before operational options ( for readablility)
        if options is None:
            print('Config not specified, using Defaults. Calc_2d_mmpad have to be added manually!')
            options = {'verbose': '0', 'max_ensemble_size': '10000000', 'tlag_domain':'time', 'nmodels':'2000000','nmod_burn': '10000',
                       'nmod_skip':'1', 'step_size': '0.05', 'term_mode':'convergence', 'conv_limit': '0.0005', 'nmod_stable': '1000',
                       'save_ensemble': '1', 'save_expvals': '1', 'nbin_1d_mppd': '50', 'nbin_2d_mppd': '50' }
        for opt in options:
            conf_SE = ElementTree.SubElement(config,opt)
            conf_SE.text=options[opt]
            if opt == 'tlag_domain':
                config.append(ElementTree.Comment(' MH Control Options '))
            elif opt == 'stepsize':
                config.append(ElementTree.Comment(' Termination Options '))
            elif opt == 'nmod_stable':
                config.append(ElementTree.Comment(' Output and Storage Options '))

        # tree.write('MTSConfig.xml')
        self._write_pretty_xml(root,file='{}/MTSConfig.xml'.format(self.opath))

    def _write_pretty_xml(self,root,file):
        """Return a pretty-printed XML string for the Element.
        """
        rough_string = ElementTree.tostring(root)
        reparsed = minidom.parseString(rough_string)
        print(reparsed)
        with open('{}'.format(file),'w') as writer :
            reparsed.writexml(writer,indent="    ",addindent="    ",newl="\n",encoding="utf-8")

        print('XML written to {}'.format(file))

    def _write_domain_counts(self,file):
        """
        Adds the domains phase counts to the domains DF and then writes it out
        """
        self.doms['Lower_Count'] = self.doms['BIN'].map(self.ldom_c)
        self.doms['Upper_Count'] = self.doms['BIN'].map(self.udom_c)
        self.doms = self.doms.fillna(0)
        self.doms.to_csv(file,sep=' ',index=False)

def phases2sdb():
    '''
    Function to take data from multiple phases and make in into a single file to be used in Pathsetting
    '''
    date_time_convert = {'TIME': lambda x: str(x),'DATE': lambda x : str(x)}
    df = pd.read_csv('E_pacific_05.pairs',converters=date_time_convert,delim_whitespace=True)
    df_snks = pd.read_csv('E_pacific_05_binned.pairs',converters=date_time_convert,delim_whitespace=True)
    cols = ['TLAG_SKS','DTLAG_SKS','FAST_SKS','DFAST_SKS','TLAG_SKKS','DTLAG_SKKS','FAST_SKKS',
            'DFAST_SKKS','SNR_SKS','SNR_SKKS']
    df_snks[cols] = df[cols]
    df_scs = pd.read_csv('Jacks_Final_ScS_overlap.rdb',converters=date_time_convert,delim_whitespace=True)
    with open('E_pacific_all_data_binned.sdb','w+') as writer:
        writer.write('STAT DATE TIME PHASE EVLA EVLO EVDP STLA STLO DIST BAZ AZI LOWMM_LAT LOWMM_LON FAST DFAST TLAG DTLAG Q SNR LOWER_BIN RSIDE_BIN \n')
        for ph in ['SKS','SKKS']:
            pplon = '{}_PP_LON'.format(ph)
            pplat = '{}_PP_LAT'.format(ph)
            q = 'Q_{}'.format(ph)
            snr = 'SNR_{}'.format(ph)
            bn = '{}_BIN'.format(ph)
            fast = 'FAST_{}'.format(ph)
            dfast = 'DFAST_{}'.format(ph)
            tlag = 'TLAG_{}'.format(ph)
            dtlag = 'DTLAG_{}'.format(ph)
            for i,row in df_snks.iterrows():
                writer.write('{} {} {} {} {:4.3f} {:4.3f} {:4.3f} {:4.3f} {:4.3f} {:4.3f} {:4.3f} {:4.3f} {:4.3f} {:4.3f} {:4.3f} {:4.3f} {:4.3f} {:4.3f} {:4.3f} {:4.3f} {} {} \n'.format(
                row['STAT'], row['DATE'], row['TIME'], ph,  row['EVLA'], row['EVLO'], row['EVDP'],
                row['STLA'],row['STLO'], row['DIST'], row['BAZ'], row['AZI'], row[pplat], row[pplon],
                row[fast], row[dfast], row[tlag], row[dtlag], row[q], row[snr], row[bn], row['RSIDE_BIN']))
        for i, row in df_scs.iterrows():
            writer.write('{} {} {} {} {:4.3f} {:4.3f} {:4.3f} {:4.3f} {:4.3f} {:4.3f} {:4.3f} {:4.3f} {:4.3f} {:4.3f} {:4.3f} {:4.3f} {:4.3f} {:4.3f} {:4.3f} {:4.3f} {} {} \n'.format(
            row['STAT'], row['DATE'], row['TIME'], 'ScS',  row['EVLA'], row['EVLO'], row['EVDP'],
            row['STLA'],row['STLO'], row['DIST'], row['BAZ'], row['AZI'], row['BNCLAT'], row['BNCLON'],
            row['FAST'], row['DFAST'], row['TLAG'], row['DTLAG'], row['Q'], row['SNR'], row['LOWMM_BIN'], row['RSIDE_BIN']))
    print('Done')

def find_phases_in_domain(phasefile,dom_ID,layer='Lower',crit=5.0,save=False):
    '''
    This function searches the input phase file for all paths that go through the given domain (intended use is D'' domains). This subset of phases is then written out for future use.
    '''
    date_time_convert = {'TIME': lambda x: str(x),'DATE': lambda x : str(x)}
    df = pd.read_csv(phasefile,converters=date_time_convert,delim_whitespace=True)
    domains = np.loadtxt('/Users/ja17375/SWSTomo/Inversions/T3_global.bins',skiprows=1)
#     dom = domains[domains[:,0] == dom_ID]
    dom = domains[np.isin(domains[:,0],dom_ID)]
    dlat = dom[0,1]
    dlon = dom[0,2]
    idx = [ ] # empty list to hold indexs of phases we want
    for i, phase in df.iterrows():
        if layer == 'Lower':
            plat = phase.LOWMM_LAT
            plon = phase.LOWMM_LON
        elif phase == 'RSide':
            plat = phase.STLA
            plon = phase.STLO

        dist = vincenty_dist(dlat, dlon, plat, plon)[0]

        if dist <= crit:
            #print(phase.STAT, phase.PHASE, dist)
            #Is the phase within a fixed distance from domain center
            idx.append(i)
    df_out = df.iloc[idx]
    if save:
        df_out.to_csv('Domain_{}_SNR10_goodQ_allphases.sdb'.format(dom_ID),index=False,sep=' ')
    else:
        return df_out

def split_by_stats(phasefile, stats):
    '''
    function to breakup phase list by station and write out these subsections

    Args:
        phasefile (str) - path to phase file to use
        stats (array-like) - array of stations to iterate over

    Returns:

        writes out textfiles of phases lists for each station
    '''
    date_time_convert = {'TIME': lambda x: str(x),'DATE': lambda x : str(x)}
    df = pd.read_csv(phasefile,converters=date_time_convert,delim_whitespace=True)

    for stat in stats:
        stat_df = df[df.STAT == stat]
        print(stat, len(stat_df))
        stat_df.to_csv('./StatCorr/{}_goodQ.sdb'.format(stat),sep=' ',index=False)

    print('Done')

if __name__ == '__main__':
  # Set = PathSetter(phasefile='/Users/ja17375/SWSTomo/Inversions/HQ_phases_on_fast_anom.sdb',
    # odir='/Users/ja17375/SWSTomo/Inversions/EP_fast_anom/')
    # Set.gen_Model_XML(stations='all', mod_name='EP_fast_anom_Model', Low_Domains=['fast_anom'])
    # Make Pathsets for Epac Data
    # Set = PathSetter(phasefile='/Users/ja17375/SWSTomo/Inversions/HQ_phases_on_fast_anom.sdb',
    # odir='/Users/ja17375/SWSTomo/Inversions/', model='EP_fast_anom_Model.xml')
    # Set.gen_PathSet_XML(stations='all',lower_dom_no='fast_anom',
    #                     phases=['ScS','SKS','SKKS'],fname='EP_fast_anom_Paths')
    # Make Pathsets for synthetics
    SynSet = PathSetter(
        phasefile='/Users/ja17375/SWSTomo/Inversions/SynthTests/RealBAZ/ppv1_real_noise10_SYNTH_sheba_results.sdb',
        odir='/Users/ja17375/SWSTomo/Inversions/SynthTests/RealBAZ/Noise10', model='Synth_Model.xml')
    SynSet.gen_synth_PathSet(lower_dom_no='fast_anom', syntest='RealBAZ/Noise10',fname='Synth_Pathset')

max_stars_repo_name: Jasplet/SWSTomo

max_issues_repo_name: Jasplet/SWSTomo

max_forks_repo_name: Jasplet/SWSTomo


text:
import matplotlib.pyplot as plt
import yaml
import os

import zipfile
import numpy as np
import torch
from torchvision import utils

# ['all_grad_norm', 'all_updated_grad_norm', 'all_probability', 'all_entropy', 'all_entropy_scale', 'all_probability_distribution']

fontdict = {
    'fontsize': 20,
    'family': 'Times New Roman',
    'weight': 'light',
}

images = {}
probability = {}
grad_norm = {}

selected_image_ids = {
    628: [19, 6, 13, 8],
    417: [0, 4, 9, 16],
    449: [2,6,23,12],
    508: [3,7,0,17]
}

os.makedirs('/workspace/mnt/storage/guangcongzheng/zju_zgc/guided-diffusion/imgs/compared_img_steps250', exist_ok=True)
save_path = '/workspace/mnt/storage/guangcongzheng/zju_zgc/guided-diffusion/imgs/compared_img_steps250/compared_img.pdf'

total_images = {
    'baseline': [],
    'ECT+EDS': []
}
for key, value in selected_image_ids.items():

    for method in ['baseline', 'ECT+EDS']:
        npz_path = '/workspace/mnt/storage/guangcongzheng/zju_zgc/guided-diffusion/log/get_compared_image/predict/model500000_steps250_scale10.0_sample50000/class{}/metainfo_{}.npz'.format(key,
                                                                                method)
        npz = np.load(npz_path)
        images[method] = npz['batch_image']  # (B,H,W,C)
        print(images[method].shape)
        for id in value:
            total_images[method].append(torch.Tensor(images[method][id]))

i = 0
j = 0

final_images = []
save_path = '/workspace/mnt/storage/guangcongzheng/zju_zgc/guided-diffusion/imgs/compared_img_steps250/{}.pdf'.format( 'baseline')
for key, value in selected_image_ids.items():

    for _ in range(4):
        print(i)
        final_images.append(total_images['baseline'][i])
        i += 1

utils.save_image(
    final_images,
    save_path,
    nrow=4,
    normalize=True,
    range=(-1, 1),
)

final_images = []
save_path = '/workspace/mnt/storage/guangcongzheng/zju_zgc/guided-diffusion/imgs/compared_img_steps250/{}.pdf'.format( 'ECT+EDS')
for key, value in selected_image_ids.items():

    for _ in range(4):
        final_images.append(total_images['ECT+EDS'][j])
        j += 1

utils.save_image(
    final_images,
    save_path,
    nrow=4,
    normalize=True,
    range=(-1, 1),
)

ext: py

max_stars_repo_name: ZGCTroy/guided-diffusion

max_issues_repo_name: ZGCTroy/guided-diffusion

max_forks_repo_name: ZGCTroy/guided-diffusion

reason: import numpy


text:
""" This module computes alignment solutions between all "a priori" solutions
    for a dataset and GAIA.

"""
import pytest

import numpy as np

from drizzlepac.haputils import testutils

from ..resources import BaseACS, BaseWFC3

def compare_apriori(dataset):
    """This test will perform fits between ALL a priori solutions and GAIA.

    Success criteria:
      * Successful fit to GAIA for dataset
      * Fit was not compromised (fit_qual != 5)
      * radial offset for new WCS < radial offset for IDC_* WCS
      * scale for new WCS within 0.1% of scale for IDC_* WCS
      * rotation for new WCS is within 0.1% of rotation from IDC_* WCS

    ASSUMPTIONS:
      * OPUS-based WCS solutions are ignored
      * Some WCS's may be defined for one image in an ASN but not another,
        in which case, that WCS is ignored (silently).

    """
    # Perform alignment of all WCS solutions with GAIA
    results_dict = testutils.compare_wcs_alignment(dataset)
    limit = 0.001

    # Now compare results to see whether the a priori solutions succeeded
    # in improving the astrometry compared to default telescope pointing
    # reported by pipeline-defined WCS (IDC_* solution)
    wcsnames = list(results_dict.keys())
    for idc_name in wcsnames:
        if 'IDC' in idc_name and '-' not in idc_name:
            wcsnames.remove(idc_name)
            break
    else:
        raise ValueError
    # Define pipeline-default fit
    pipeline_results = results_dict[idc_name]
    pipeline_offset = np.sqrt(pipeline_results['offset_x']**2 +
                              pipeline_results['offset_y']**2)

    success = False
    # compare each WCS fit results with pipeline-default fit
    for wcs in wcsnames:
        print("Comparing fit for WCS='{}'...".format(wcs), end=' ')
        results = results_dict[wcs]

        # Check that fit for this WCS was successful
        status = (results['status'] == 0).all()
        # check that fit was not compromised or otherwise invalid
        fit_qual = (results['fit_qual'] < 5).all()

        # Check radial offset for this WCS compared to radial offset for
        # IDC* WCS
        offset = np.sqrt(results['offset_x']**2 + results['offset_y']**2)
        delta = ((offset < pipeline_offset).all() or
                 np.allclose(offset, pipeline_offset, rtol=0, atol=1))

        # Check that rotation and scale are within
        rot = np.allclose(results['rotation'], pipeline_results['rotation'],
                          rtol=limit, atol=0)
        scale = np.allclose(results['scale'], pipeline_results['scale'],
                            rtol=limit, atol=0)

        # Determine success/failure of this dataset's fit
        if all([status, fit_qual, delta, rot, scale]):
            # If at least one WCS succeeds, overall success
            # needs to be set to True
            success = True
            print("SUCCESS")
        else:
            print("FAILED  due to:")
            if not status:
                print("\t* invalid STATUS")
            if not fit_qual:
                print("\t* invalid fit quality")
            if not delta:
                print("\t* increased offset from GAIA.")
            if not rot:
                print("\t* larger rotation from fit.")
            if not scale:
                print("\t* larger scale from fit.")

    assert success

class TestAcsApriori(BaseACS):
    """ Tests which validate whether mosaics can be aligned to an astrometric standard,
        evaluate the quality of the fit, and generate a new WCS.

        * These can only be run using the pytest option:

            pytest --bigdata

        * This test is also compatible with pytest-xdist.
    """

    @pytest.mark.bigdata
    @pytest.mark.parametrize('dataset', ['jb1601020', 'J9I408010'])
    def test_apriori(self, dataset):
        compare_apriori(dataset)

class TestWFC3Apriori(BaseWFC3):
    """ Tests which validate whether mosaics can be aligned to an astrometric
        standard, evaluate the quality of the fit, and generate a new WCS.

        * These can only be run using the pytest option:

            pytest --bigdata

        * This test is also compatible with pytest-xdist.

    """

    @pytest.mark.bigdata
    @pytest.mark.parametrize(
        'dataset', ['ic0g0l010', 'icnw34040']
    )
    def test_apriori(self, dataset):
        compare_apriori(dataset)

max_stars_repo_name: check-spelling/drizzlepac

max_issues_repo_name: check-spelling/drizzlepac

max_forks_repo_name: check-spelling/drizzlepac


text:
from abc import ABC
from functools import lru_cache
from inspect import getfullargspec
from itertools import chain, repeat
from typing import \
    (
    Tuple, Dict, Union, Set, Callable, List
)

import numpy as np

from autofit.graphical.factor_graphs.abstract import \
    (
    AbstractNode, FactorValue
)
from autofit.graphical.utils import \
    (
    aggregate, Axis, cached_property
)
from autofit.mapper.variable import Variable

class AbstractFactor(AbstractNode, ABC):
    def __init__(
            self,
            name="",
            **kwargs: Variable,
    ):
        super().__init__(**kwargs)
        self._name = name or f"factor_{self.id}"
        self._deterministic_variables = set()

    @property
    def deterministic_variables(self) -> Set[Variable]:
        return self._deterministic_variables

    @property
    def name(self) -> str:
        return self._name

    def __mul__(self, other):
        """
        When two factors are multiplied together this creates a graph
        """
        from autofit.graphical.factor_graphs.graph import FactorGraph
        return FactorGraph([self]) * other

    @property
    def variables(self) -> Set[Variable]:
        """
        Dictionary mapping the names of variables to those variables
        """
        return set(self._kwargs.values())

    @property
    def _kwargs_dims(self) -> Dict[str, int]:
        """
        The number of plates for each keyword argument variable
        """
        return {
            key: len(value)
            for key, value
            in self._kwargs.items()
        }

    @cached_property
    def _variable_plates(self) -> Dict[str, np.ndarray]:
        """
        Maps the name of each variable to the indices of its plates
        within this node
        """
        return {
            variable: self._match_plates(
                variable.plates
            )
            for variable
            in self.all_variables
        }

    @property
    def n_deterministic(self) -> int:
        """
        How many deterministic variables are there associated with this node?
        """
        return len(self._deterministic_variables)

    def __hash__(self):
        return hash((type(self), self.id))

    def _resolve_args(
            self,
            **kwargs: np.ndarray
    ) -> dict:
        """
        Transforms in the input arguments to match the arguments
        specified for the factor.

        Parameters
        ----------
        args
        kwargs

        Returns
        -------

        """
        return {n: kwargs[v.name] for n, v in self._kwargs.items()}

class Factor(AbstractFactor):
    """
    A node in a graph representing a factor with analytic evaluation
    of its Jacobian

    Parameters
    ----------
    factor
        the function being wrapped, must accept calls through keyword
        argument

    name: optional, str
        the name of the factor, if not passed then uses the name
        of the function passed

    vectorised: optional, bool
        if true the factor will call the function directly over multiple
        inputs. If false the factor will call the function iteratively
        over each argument.

    is_scalar: optional, bool
        if true the factor returns a scalar value. Note if multiple arguments
        are passed then a vector will still be returned

    kwargs: Variables
        Variables for each keyword argument for the function

    Methods
    -------
    __call__({x: x0}, axis=axis)
        calls the factor, the passed input must be a dictionary with
        where the keys are the Variable objects that the function takes
        as input. The Variable keys only have to match the _names_
        of the variables of the function.

        `axis` controls the shape of the output if the variables and factor have
        plates associated with them, when axis=False then no reduction is performed,
        otherwise it is equivalent to calling np.sum(log_val, axis=axis) on the
        returned value

        returns a FactorValue object which behaves like an np.ndarray

    func_jacobian({x: x0}, variables=(x,), axis=axis)
        calls the factor and returns it value and the jacobian of its value
        with respect to the `variables` passed. if variables is None then
        it returns the jacobian with respect to all variables.

        returns fval, {x: d fval / dx}
    """

    def __init__(
            self,
            factor: Callable,
            name="",
            vectorised=False,
            is_scalar=False,
            **kwargs: Variable
    ):
        """
        A node in a graph representing a factor

        Parameters
        ----------
        factor
            A wrapper around some callable
        args
            Variables representing positional arguments for the function
        kwargs
            Variables representing keyword arguments for the function
        """
        self.vectorised = vectorised
        self.is_scalar = is_scalar
        self._set_factor(factor)

        args = getfullargspec(self._factor).args
        kwargs = {
            **kwargs,
            **{
                arg: Variable(arg)
                for arg
                in args
                if arg not in kwargs and arg != "self"
            }
        }

        super().__init__(
            **kwargs,
            name=name or factor.__name__
        )

    def _set_factor(self, factor):
        self._factor = factor
        self._has_exact_projection = getattr(
            factor, 'has_exact_projection', None)
        self._calc_exact_projection = getattr(
            factor, 'calc_exact_projection', None)
        self._calc_exact_update = getattr(
            factor, 'calc_exact_update', None)

    def has_exact_projection(self, mean_field) -> bool:
        if self._has_exact_projection:
            return self._has_exact_projection(
                **self.resolve_variable_dict(mean_field))
        else:
            return False

    def calc_exact_projection(self, mean_field) -> 'MeanField':
        if self._calc_exact_projection:
            from autofit.graphical.mean_field import MeanField
            projection = self._calc_exact_projection(
                **self.resolve_variable_dict(mean_field))
            return MeanField({
                self._kwargs[v]: dist for v, dist in projection.items()
            })
        else:
            return NotImplementedError

    def calc_exact_update(self, mean_field) -> 'MeanField':
        if self._calc_exact_update:
            from autofit.graphical.mean_field import MeanField
            projection = self._calc_exact_update(
                **self.resolve_variable_dict(mean_field))
            return MeanField({
                self._kwargs[v]: dist for v, dist in projection.items()
            })
        else:
            return NotImplementedError

    def safe_exact_update(self, mean_field) -> Tuple[bool, 'MeanField']:
        if self._has_exact_projection:
            from autofit.graphical.mean_field import MeanField
            _mean_field = self.resolve_variable_dict(mean_field)
            if self._has_exact_projection(**_mean_field):
                projection = self._calc_exact_update(**_mean_field)
                return True, MeanField({
                    self._kwargs[v]: dist for v, dist in projection.items()
                })

        return False, mean_field

    def __hash__(self) -> int:
        # TODO: might this break factor repetition somewhere?
        return hash(self._factor)

    def _reshape_factor(
            self, factor_val, values
    ):
        shift, shape = self._function_shape(**values)
        if self.is_scalar:
            if shift:
                return np.sum(
                    factor_val, axis=np.arange(1, np.ndim(factor_val)))
            return np.sum(factor_val)
        return np.reshape(factor_val, shape)

    def _function_shape(
            self,
            **kwargs: np.ndarray) -> Tuple[int, ...]:
        """
        Calculates the expected function shape based on the variables
        """
        var_shapes = {
            k: np.shape(x) for k, x in kwargs.items()}
        return self._var_shape(**var_shapes)

    @lru_cache()
    def _var_shape(self, **kwargs: Tuple[int, ...]) -> Tuple[int, ...]:
        """This is called by _function_shape

        caches result so that does not have to be recalculated each call

        lru_cache caches f(x=1, y=2) to f(y=2, x=1), but in this case
        it should be find as the order of kwargs is set by self._kwargs
        which should be stable
        """
        var_shapes = {self._kwargs[k]: v for k, v in kwargs.items()}
        var_dims_diffs = {
            v: len(s) - v.ndim
            for v, s in var_shapes.items()
        }
        """
        If all the passed variables have an extra dimension then
        we assume we're evaluating multiple instances of the function at the
        same time

        otherwise an error is raised
        """
        if set(var_dims_diffs.values()) == {1}:
            # Check if we're passing multiple values e.g. for sampling
            shift = 1
        elif set(var_dims_diffs.values()) == {0}:
            shift = 0
        else:
            raise ValueError("dimensions of passed inputs do not match")

        """
        Updating shape of output array to match input arrays

        singleton dimensions are always assumed to match as in
        standard array broadcasting

        e.g. (1, 2, 3) == (3, 2, 1)
        """
        shape = np.ones(self.ndim + shift, dtype=int)
        for v, vs in var_shapes.items():
            ind = self._variable_plates[v] + shift
            vshape = vs[shift:]
            if shift:
                ind = np.r_[0, ind]
                vshape = (vs[0],) + vshape

            if shape.size:
                if not (
                        np.equal(shape[ind], 1) |
                        np.equal(shape[ind], vshape) |
                        np.equal(vshape, 1)).all():
                    raise AssertionError(
                        "Shapes do not match"
                    )
                shape[ind] = np.maximum(shape[ind], vshape)

        return shift, tuple(shape)

    def _call_factor(
            self,
            **kwargs: np.ndarray
    ) -> np.ndarray:
        """
        Call the underlying function

        Parameters
        ----------
        args
            Positional arguments for the function
        kwargs
            Keyword arguments for the function

        Returns
        -------
        Value returned by the factor
        """
        # kws = self._resolve_args(
        #     **kwargs
        # )

        if self.vectorised:
            return self._factor(**kwargs)

        """Some factors may not be vectorised to broadcast over
        multiple inputs

        this method checks whether multiple input values have been
        passed, and if so automatically loops over the inputs.
        If any of the inputs have initial dimension one, it repeats
        that value to match the length of the other inputs

        If the other inputs do not match then it raises ValueError
        """
        kwargs_dims = {k: np.ndim(a) for k, a in kwargs.items()}
        # Check dimensions of inputs directly match plates
        direct_call = (
            all(dim == kwargs_dims[k] for k, dim in self._kwargs_dims.items()))
        if direct_call:
            return self._factor(**kwargs)

        # Check dimensions of inputs match plates + 1
        vectorised = (
            all(dim + 1 == kwargs_dims[k]
                for k, dim in self._kwargs_dims.items()))

        if not vectorised:
            raise ValueError(
                "input dimensions do not match required dims"
                f"input: **kwargs={kwargs_dims}"
                f"required: "
                f"**kwargs={self._kwargs_dims}")

        kw_lens = {k: len(a) for k, a in kwargs.items()}

        # checking 1st dimensions match
        sizes = set(kw_lens.values())
        dim0 = max(sizes)
        if sizes.difference({1, dim0}):
            raise ValueError(
                f"size mismatch first dimensions passed: {sizes}")

        iter_kws = {
            k: iter(a) if kw_lens[k] == dim0 else iter(repeat(a[0]))
            for k, a in kwargs.items()}

        # iterator to generate keyword arguments
        def gen_kwargs():
            for _ in range(dim0):
                yield {
                    k: next(a) for k, a in iter_kws.items()}

        # TODO this loop can also be parallelised for increased performance
        res = np.array([
            self._factor(**kws)
            for kws in gen_kwargs()])

        return res

    def __call__(
            self,
            variable_dict: Dict[
                Variable,
                Union[
                    np.ndarray,
                    float,
                    List[float]
                ]
            ],
            axis: Axis = False,
    ) -> FactorValue:
        """
        Call the underlying factor

        Parameters
        ----------
        args
            Positional arguments for the factor
        kwargs
            Keyword arguments for the factor

        Returns
        -------
        Object encapsulating the result of the function call
        """
        kwargs = self.resolve_variable_dict(variable_dict)
        val = self._call_factor(**kwargs)
        val = aggregate(self._reshape_factor(val, kwargs), axis)
        return FactorValue(val, {})

    def broadcast_variable(
            self,
            variable: str,
            value: np.ndarray
    ) -> np.ndarray:
        """
        broadcasts the value of a variable to match the specific shape
        of the factor

        if the number of dimensions passed of the variable is 1
        greater than the dimensions of the variable then it's assumed
        that that dimension corresponds to multiple samples of that variable
        """
        return self._broadcast(
            self._variable_plates[variable],
            value
        )

    def collapse(
            self,
            variable: str,
            value: np.ndarray,
            agg_func=np.sum
    ) -> np.ndarray:
        """
        broadcasts `value` to match the specific shape of the factor,
        where `value` has the shape of the factor

        if the number of dimensions passed of the variable is 1
        greater than the dimensions of the variable then it's assumed
        that that dimension corresponds to multiple samples of that variable
        """
        ndim = np.ndim(value)
        shift = ndim - self.ndim
        assert shift in {0, 1}
        inds = self._variable_plates[variable] + shift
        dropaxes = tuple(np.setdiff1d(
            np.arange(shift, ndim), inds))

        # to ensured axes of returned array is in the correct order
        moved = np.moveaxis(value, inds, np.sort(inds))
        return agg_func(moved, axis=dropaxes)

    def __eq__(self, other: Union["Factor", Variable]):
        """
        If set equal to a variable that variable is taken to be deterministic and
        so a DeterministicFactorNode is generated.
        """
        if isinstance(other, Factor):
            if isinstance(other, type(self)):
                return (
                        (self._factor == other._factor)
                        and (frozenset(self._kwargs.items())
                             == frozenset(other._kwargs.items()))
                        and (frozenset(self.variables)
                             == frozenset(other.variables))
                        and (frozenset(self.deterministic_variables)
                             == frozenset(self.deterministic_variables)))
            else:
                return False

        return DeterministicFactor(
            self._factor,
            other,
            name=self.name,
            **self._kwargs
        )

    def __repr__(self) -> str:
        args = ", ".join(chain(
            map("{0[0]}={0[1]}".format, self._kwargs.items())))
        return f"Factor({self.name}, {args})"

class DeterministicFactor(Factor):
    """
    A deterministic factor is used to convert a function f(g(x)) to f(y)g(x) (integrating over y wit
    a delta function) so that it can be represented in a factor graph.

    Parameters
    ----------
    factor
        The original factor to which the deterministic factor is associated
    variable
        The deterministic variable that is returned by the factor, so
        to represent the case f(g(x)), we would define,

        ```
        >>> x = Variable('x')
        >>> y = Variable('y')
        >>> g_ = Factor(g, x) == y
        >>> f_ = Factor(f, y)
        ```
        Alternatively g could be directly defined,
        ```
        >>> g_ = DeterministicFactor(g, y, x=x)
        ```

    kwargs
        Variables for the original factor
    """

    def __init__(
            self,
            factor: Callable,
            variable: Variable,
            *args: Variable,
            name: str = '',
            **kwargs: Variable
    ):
        """
        A deterministic factor is used to convert a function f(g(x)) to f(y)g(x) (integrating over y wit
        a delta function) so that it can be represented in a factor graph.

        Parameters
        ----------
        factor
            The original factor to which the deterministic factor is associated
        variable
            The deterministic factor used
        args
            Variables for the original factor
        kwargs
            Variables for the original factor
        """
        super().__init__(
            factor,
            *args,
            name=name or factor.__name__,
            **kwargs
        )
        self._deterministic_variables = {
            variable
        }

    def __call__(
            self,
            variable_dict: Dict[Variable, np.ndarray],
            axis: Axis = False,
            # **kwargs: np.ndarray
    ) -> FactorValue:
        """
        Call this factor with a set of arguments

        Parameters
        ----------
        args
            Positional arguments for the underlying factor
        kwargs
            Keyword arguments for the underlying factor

        Returns
        -------
        An object encapsulating the value for the factor
        """
        kwargs = self.resolve_variable_dict(variable_dict)
        res = self._call_factor(**kwargs)
        shift, shape = self._function_shape(**kwargs)
        plate_dim = dict(zip(self.plates, shape[shift:]))

        det_shapes = {
            v: shape[:shift] + tuple(
                plate_dim[p] for p in v.plates)
            for v in self.deterministic_variables
        }

        if not (isinstance(res, tuple) or self.n_deterministic > 1):
            res = res,

        log_val = (
            0. if (shape == () or axis is None) else
            aggregate(np.zeros(tuple(1 for _ in shape)), axis))
        det_vals = {
            k: np.reshape(val, det_shapes[k])
            if det_shapes[k]
            else val
            for k, val
            in zip(self._deterministic_variables, res)
        }
        return FactorValue(log_val, det_vals)

    def __repr__(self) -> str:
        factor_str = super().__repr__()
        var_str = ", ".join(sorted(variable.name for variable in self._deterministic_variables))
        return f"({factor_str} == ({var_str}))"

max_stars_repo_name: vishalbelsare/PyAutoFit

max_issues_repo_name: vishalbelsare/PyAutoFit

max_forks_repo_name: vishalbelsare/PyAutoFit


text:
import numpy as np
import tensorflow as tf
from tensorflow import keras
from tensorflow.keras.applications import mobilenet_v2
from tensorflow.python.keras.preprocessing.image import ImageDataGenerator
from tensorflow.python.training.basic_session_run_hooks import CheckpointSaverListener

class MyListener(CheckpointSaverListener):
    def begin(self):
        print(".")

    def end(self, session, global_step_value):
        print("_")

if __name__ == '__main__':
    base_path = "../../../data"
    output_path = "../../../output"
    images_dir = "%s/flower_photos" % base_path

    #####################
    # Data Generator
    #####################

    # Before do this, please add new classes at the end of ImageNetLabels.txt. each class per line
    # classes: daisy、dandelion、roses、sunflowers、tulips
    labels_path = "ImageNetLabels.txt"
    imagenet_labels = np.array(open(labels_path).read().splitlines())[1:].tolist()  # drop first line 'background'

    print(len(imagenet_labels), imagenet_labels)

    train_datagen = ImageDataGenerator(preprocessing_function=mobilenet_v2.preprocess_input,
                                       validation_split=0.2)
    train_generator = train_datagen.flow_from_directory(images_dir, target_size=(224, 224),
                                                        color_mode='rgb', batch_size=32, shuffle=True,
                                                        classes=imagenet_labels, class_mode='sparse',
                                                        subset='training')
    validation_generator = train_datagen.flow_from_directory(images_dir, target_size=(224, 224),
                                                             color_mode='rgb', batch_size=32, shuffle=True,
                                                             classes=imagenet_labels, class_mode='sparse',
                                                             subset='validation')

    ###################
    # Define transfer learning model
    ###################

    base_model = mobilenet_v2.MobileNetV2(weights='imagenet', include_top=False, input_shape=(224, 224, 3))
    base_model.trainable = False

    keras_model = keras.Sequential([
        base_model,
        keras.layers.GlobalAveragePooling2D(),
        keras.layers.Dense(1024, activation='relu'),
        keras.layers.Dense(512, activation='relu'),
        keras.layers.Dropout(0.2),
        keras.layers.Dense(5, activation='softmax')
    ])
    keras_model.compile(optimizer='adam', loss='sparse_categorical_crossentropy', metric='accuracy')

    ###################
    # Use tensorflow estimator
    ###################

    # Multi GPUs
    strategy = tf.contrib.distribute.MirroredStrategy()
    config = tf.estimator.RunConfig(train_distribute=strategy)

    est_model = keras.estimator.model_to_estimator(keras_model=keras_model, config=config)

    train_x = []
    train_y = []
    for inputs_batch, labels_batch in validation_generator:
        train_x.append(inputs_batch)
        train_y.append(labels_batch)

    train_input_fn = tf.estimator.inputs.numpy_input_fn(x=train_x, y=train_y, num_epochs=1, shuffle=False)

    my_listener = MyListener()

    # To train, we call Estimator's train function:
    est_model.train(input_fn=train_input_fn, steps=5, saving_listeners=my_listener)


text:
#!/usr/bin/env python3
# Original implementation by Federico Magistri
import os

import cv2 as cv
import numpy as np
import open3d as o3d

class ScanNet:
    def __init__(self, data_source, get_color=False, split=None):

        self.name = data_source.split("/")[-1]

        self.get_color = get_color
        self.data_source = data_source

        self.rgb_list = os.listdir(os.path.join(self.data_source, "color"))
        self.rgb_list.sort()

        self.depth_list = os.listdir(os.path.join(self.data_source, "depth"))
        self.depth_list.sort()

        self.pose_list = os.listdir(os.path.join(self.data_source, "pose"))
        self.pose_list.sort()

    def __len__(self):
        return len(self.pose_list)

    def __getitem__(self, idx):

        pose = np.loadtxt(os.path.join(self.data_source, "pose", self.pose_list[idx]))
        pose_inv = np.linalg.inv(pose)

        bgr = cv.imread(os.path.join(self.data_source, "color", self.rgb_list[idx]))
        rgb = cv.cvtColor(bgr, cv.COLOR_BGR2RGB)

        depth = cv.imread(os.path.join(self.data_source, "depth", self.depth_list[idx]), -1)

        K = np.loadtxt(os.path.join(self.data_source, "intrinsic_depth.txt"))

        rgb = o3d.geometry.Image(rgb)
        depth = o3d.geometry.Image(depth)

        h, w, _ = bgr.shape

        intrinsic = o3d.camera.PinholeCameraIntrinsic()
        intrinsic.set_intrinsics(
            height=h,
            width=w,
            fx=K[0, 0],
            fy=K[1, 1],
            cx=K[0, 2],
            cy=K[1, 2],
        )

        rgbd = o3d.geometry.RGBDImage.create_from_color_and_depth(
            rgb, depth, depth_scale=1000, depth_trunc=1000.0, convert_rgb_to_intensity=False
        )

        pcd = o3d.geometry.PointCloud.create_from_rgbd_image(
            rgbd, intrinsic, pose_inv, project_valid_depth_only=True
        )

        xyz = np.array(pcd.points)
        colors = np.array(pcd.colors)

        xyz = np.array(xyz)
        colors = np.array(colors)

        if self.get_color:
            return xyz, colors, np.array(pose)

        return xyz, np.array(pose)

max_stars_repo_name: Willyzw/vdbfusion

max_issues_repo_name: Willyzw/vdbfusion

max_forks_repo_name: Willyzw/vdbfusion


text:
import numpy as np
import random
import os

import torch

from torch.utils.data.sampler import Sampler
from torch import optim
from torch.autograd import Variable

from configs import g_conf

# TODO: When putting sequences, the steering continuity and integrity needs to be verified
def get_rank(input_array):

    rank = 0
    while True:
        try:
            length = len(input_array)
            input_array = input_array[0]
            rank += 1
        except:
            return rank

class RandomSampler(Sampler):
    r"""Samples elements randomly from a given list

    Arguments:
        indices (list): a list of indices
    """

    def __init__(self, keys, executed_iterations):
        self.iterations_to_execute = ((g_conf.NUMBER_ITERATIONS) * g_conf.BATCH_SIZE) -\
                                     (executed_iterations)

        self.keys = keys

    def __iter__(self):

        return iter([random.choice(self.keys) for _ in range(self.iterations_to_execute)])

    def __len__(self):
        return self.iterations_to_execute

class SubsetSampler(Sampler):
    r"""Samples elements randomly from a given list of indices, without replacement.

    Arguments:
        indices (list): a list of indices
    """

    def __init__(self, indices):
        self.indices = indices

    def __iter__(self):
        return (self.indices[i] for i in range(len(self.indices)))

    def __len__(self):
        return len(self.indices)

class PreSplittedSampler(Sampler):
    """ Sample on a list of keys that was previously splitted

    """

    def __init__(self, keys, executed_iterations, weights=None):

        self.keys = keys
        if (weights is None) or (weights == 'None'):
            self.weights = np.asarray([1.0/float(len(self.keys))]*len(self.keys), dtype=np.float)
        else:
            self.weights = np.asarray(weights)

        self.iterations_to_execute = g_conf.NUMBER_ITERATIONS * g_conf.BATCH_SIZE -\
                                     executed_iterations + g_conf.BATCH_SIZE
        self.replacement = True

    def __iter__(self):
        """

            OBS: One possible thing to be done is the possibility to have a matrix of ids
            of rank N
            OBS2: Probably we dont need weights right now

        Returns:
            Iterator to get ids for the dataset

        """

        rank_keys = get_rank(self.keys)

        # First we check how many subdivisions there are
        weights = torch.from_numpy(self.weights)

        if rank_keys == 2:
            idx = torch.multinomial(weights, self.iterations_to_execute, True)
            idx = idx.tolist()
            return iter([random.choice(self.keys[i]) for i in idx])

        elif rank_keys == 3:
            weights = torch.tensor([1.0 / float(len(self.keys))] * len(self.keys),
                                   dtype=torch.double)
            idx = torch.multinomial(weights, self.iterations_to_execute, True)
            idx = idx.tolist()
            weights = torch.tensor([1.0 / float(len(self.keys[0]))] * len(self.keys[0]),
                                   dtype=torch.double)
            idy = torch.multinomial(weights, self.iterations_to_execute, True)
            idy = idy.tolist()

            return iter([random.choice(self.keys[i][j]) for i, j in zip(idx,idy)])

        else:
            raise ValueError("Keys have invalid rank")

    def __len__(self):
        return self.iterations_to_execute

class LogitSplittedSampler(Sampler):
    """ Sample on a list of keys that was previously splitted
        weights for sampling are logits and have to be softmaxed

    """

    def __init__(self, keys, executed_iterations, weights=None):

        self.keys = keys
        if weights is None:
            self.weights = torch.tensor([1.0/float(len(self.keys))]*len(self.keys), dtype=torch.double)
        else:
            self.weights = torch.from_numpy(weights)
        self.weights = Variable(self.weights, requires_grad=True)

        assert len(self.weights) == len(self.keys), "Number of weights and keys should be the same"
        self.iterations_to_execute = g_conf.NUMBER_ITERATIONS * g_conf.BATCH_SIZE -\
                                     executed_iterations + g_conf.BATCH_SIZE
        self.replacement = True
        self.optim = optim.Adam([self.weights,], lr=0.01)

    def __iter__(self):
        """

            OBS: One possible thing to be done is the possibility to have a matrix of ids
            of rank N
            OBS2: Probably we dont need weights right now

        Returns:
            Iterator to get ids for the dataset

        """
        weights = F.softmax(self.weights)
        idx = torch.multinomial(weights, self.iterations_to_execute, True)
        idx = idx.tolist()
        return iter([random.choice(self.keys[i]) for i in idx])

    def update_weights(self, advantage, perturb=False):
        self.optim.zero_grad()
        obj = torch.sum(-self.weights * advantage)
        obj.backward()
        self.optim.step()
        if perturb:
            N = len(self.weights)
            self.weights = self.weights + torch.normal(torch.zeros(N), perturb * torch.ones(N)).double()

    def __len__(self):
        return self.iterations_to_execute

class BatchSequenceSampler(object):
    r"""Wraps another sampler to yield a mini-batch of indices taking a certain sequence size

    Args:
        sampler (Sampler): Base sampler.
        batch_size (int): Size of mini-batch.
        drop_last (bool): If ``True``, the sampler will drop the last batch if
            its size would be less than ``batch_size``

    """

    def __init__(self, keys, executed_iterations,
                 batch_size, sequence_size, sequence_stride, drop_last=True):
        sampler = PreSplittedSampler(keys, executed_iterations)

        if not isinstance(sampler, Sampler):
            raise ValueError("sampler should be an instance of "
                             "torch.utils.data.Sampler, but got sampler={}"
                             .format(sampler))
        if not isinstance(sequence_size, int) or isinstance(sequence_size, bool) or \
                sequence_size <= 0:
            raise ValueError("batch_size should be a positive integeral value, "
                             "but got batch_size={}".format(sequence_size))
        if not isinstance(drop_last, bool):
            raise ValueError("drop_last should be a boolean value, but got "
                             "drop_last={}".format(drop_last))
        self.sampler = sampler
        self.sequence_size = sequence_size
        self.batch_size = batch_size
        self.drop_last = drop_last
        self.sequence_stride = sequence_stride

    def __iter__(self):

        batch = []
        for idx in self.sampler:
            for seq in range(0, self.sequence_size * self.sequence_stride, self.sequence_stride):
                batch.append(int(idx)+seq)
            if len(batch) == self.batch_size:
                yield batch
                batch = []
        if len(batch) > 0 and not self.drop_last:
            yield batch

        #if len(batch) > 0 and not self.drop_last:
        #    yield batch

    def __len__(self):
        if self.drop_last:
            return len(self.sampler) // self.sequence_size
        else:
            return (len(self.sampler) + self.sequence_size - 1) // self.sequence_size


text:
#!/usr/bin/env python

import numpy as np
import matplotlib
import matplotlib.pyplot as plt
from scipy.signal import windows
from astropy.io import fits
from astropy import wcs

arcsec = np.pi / (180.0 * 3600)  # [radians]  = 1/206265 radian/arcsec

image_raw = plt.imread("alma-logo-white-outline.png")
print("shape", image_raw.shape)
# (M, N, 4) for RGBA images.
# (703, 498, 4)
ny, nx, nchan = image_raw.shape

image_raw = image_raw[
    :, :, 0
]  # loose the last color dimension. Since it's white, it doesn't matter which channel we take.

# Since we'll be using the FFT, and it assumes a circular image, our first step is to pad the image with zeros to make the edges consistent with each other. We'll also take this opportunity to make it a square image and flip it to conform to our $l,m$ and $u,v$ conventions, since we'll be using the FFT. These are that East ($l$) and North ($m$) *increase* with array index.

# do the padding
# ncol > nrow
# pad top and bottom to make square
N = 800
x0 = (N - nx) // 2
x1 = x0 + nx
y0 = (N - ny) // 2
y1 = y0 + ny
# now pad into a 256 x 256 array
image_pad = np.zeros((N, N), dtype="float")
image_pad[y0:y1, x0:x1] = image_raw

# flip the image to be N/S
# AND
# flip image to be E/W
# l,m are increasing with array index. l corresponds to east.
image_stock = np.fliplr(np.flipud(image_pad))

# ## Some thoughts on the information content of this image
#
# The array configuration and observing frequency specify a maxmimum spatial frequency that we sample in our measurement set. This defines the "minimum recoverable image scale" that is probed by our dataset. Now, of course, if we wanted to make prior assumptions then we could go further, but let's just for a moment consider the case where we have perfect (noiseless) Nyquist sampling.

# The Nyquist frequency is
#
# $$
# u = \frac{1}{2 \times \mathrm{cell\_size}}
# $$
#
# This "stock" image contains 800x800 pixels. The next step of the image generation process will be to assign a `cell_size` to each pixel, essentially specifying how big the image will appear on sky. Since we're working with ALMA to study protoplanetary disks, we will say that this image corresponds to a diameter of 4 arcseconds.
#
# We probably also want an image that has max spatial frequencies closer to what ALMA can actually sample. What is the maximum spatial frequency probed by 16 km baselines at 1.3mm?

max_baseline_Ml = 1e-6 * 16 / (1.3e-6)  # [kl]
print(
    "max ALMA spatial freq 16km baselines at 1.3mm is {:.3f} Mlambda".format(
        max_baseline_Ml
    )
)

# Calculate the maximum spatial frequency included in a certain sized image
def get_max_spatial_freq(npix, diameter=4):
    cell_size = diameter / npix  # arcsec
    return cell_size, 1e-6 / (2 * cell_size * arcsec)  # Mlambda

cell_size, spatial_freq = get_max_spatial_freq(512)
print(
    "Target cell_size {:} arcsec delivers max spatial freq {:}".format(
        cell_size, spatial_freq
    )
)

# a 512x512 image spanning 4 arcseconds in diameter (cell_size approx 0.008 arcsec) is a nice image size providing a maximal amount of information content that could be measured by ALMA in its longest baseline configuration at Band 6.

# # Rescaling and Resampling an Image to a lower resolution
#
# Say we'd like to take our 800x800 image and rescale it to a 512x512 image. What's the right thing to do to? [Wikipedia](https://en.wikipedia.org/wiki/Image_scaling) actually has a few helpful notes about upscaling. There's also the idea about low-pass filtering. I think really what we want to do is [low-pass filter](https://en.wikipedia.org/wiki/Low-pass_filter) and resample. This is very much looking at what the response is going to be like in the Fourier domain vs image domain. An "ideal" low pass filter has the nasty Sinc effect. We want something that has better sidelobes. First, let's take a look at what the FFT looks like (amplitude only).

fft_stock = np.fft.fftshift(np.fft.fft2(np.fft.fftshift(image_stock)))
amp_stock = np.abs(fft_stock)

def create_window(window_function, npix_start, npix_end):
    assert npix_start > npix_end, "must downscale image"
    npad = npix_start - npix_end
    assert npad % 2 == 0, "must have even numbers of start and end pixels"
    w_middle = window_function(
        npix_end, False
    )  # we want the max to be at the x = 0 pixel
    w_stub = np.zeros(npad // 2)
    return np.concatenate([w_stub, w_middle, w_stub])

def truncate_array(array, npix_start, npix_end):
    npad = npix_start - npix_end
    ntrunc = npad // 2
    return array[ntrunc:-ntrunc, ntrunc:-ntrunc]

def create_apodizer(window_function, npix_start, npix_end):
    """
    Make an image
    """
    one_D = create_window(window_function, npix_start, npix_end)
    return np.outer(one_D, one_D)

apod_blackman = create_apodizer(windows.blackmanharris, 800, 512)
apod_boxcar = create_apodizer(windows.boxcar, 800, 512)

# So the answer is clearly that apodization is very important! But, these images still have 800 x 800 pixels, what happens if we drop things down to 512x512. How would we do this? My first thought is to just get rid of the extra space in the Fourier Plane, this is equivalent to "unpadding" with zeros.

fft_blackman = fft_stock * apod_blackman
fft_blackman_trunc = truncate_array(fft_blackman, 800, 512)
image_blackman = np.fft.fftshift(np.fft.ifft2(np.fft.fftshift(fft_blackman_trunc)))

amp_blackman = np.abs(fft_blackman)

# write the FITS file
# create the FITS file
w = wcs.WCS(naxis=2)

w.wcs.crpix = np.array([1, 1])
w.wcs.cdelt = (
    np.array([-cell_size, cell_size]) * arcsec * 180.0 / np.pi
)  # decimal degrees
w.wcs.ctype = ["RA---TAN", "DEC--TAN"]

header = w.to_header()

hdu = fits.PrimaryHDU(np.fliplr(np.real(image_blackman)), header=header)
hdul = fits.HDUList([hdu])
hdul.writeto("logo_cont.fits", overwrite=True)
hdul.close()

max_stars_repo_name: MPoL-dev/mpoldatasets

max_issues_repo_name: MPoL-dev/mpoldatasets

max_forks_repo_name: MPoL-dev/mpoldatasets

reason: import numpy,from scipy,from astropy


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

import unittest

import numpy as np

from tensorflowjs import quantization

class TestQuantizationUtil(unittest.TestCase):

  def assertDictContainsSubsetAlmostEqual(self, d1, d2):
    self.assertIsInstance(d1, dict)
    self.assertIsInstance(d2, dict)

    d1_keys = set(d1.keys())
    d2_keys = set(d2.keys())

    self.assertTrue(d2_keys.issubset(d1_keys))

    for key in d2_keys:
      self.assertAlmostEqual(d1[key], d2[key])

  def _runQuantizeTest(
      self, range_min, range_max, data_dtype, quantization_dtype,
      expected_metadata):
    d = np.arange(range_min, range_max + 1, dtype=data_dtype)
    q, metadata = quantization.quantize_weights(d, quantization_dtype)

    self.assertDictContainsSubsetAlmostEqual(metadata, expected_metadata)
    self.assertEqual(q.dtype, quantization_dtype)

    de_q = quantization.dequantize_weights(
        q, metadata, data_dtype)
    if data_dtype != np.float32:
      np.testing.assert_allclose(de_q, d)
    else:
      np.testing.assert_array_almost_equal(de_q, d, decimal=2)

    if quantization_dtype in [np.uint8, np.uint16]:
      s = metadata['scale']
      m = metadata['min']
      if range_min <= 0 <= range_max:
        d_0 = np.zeros(1, data_dtype)
        q_0 = np.round((d_0 - m) / s).astype(quantization_dtype)
        self.assertEqual(
            quantization.dequantize_weights(q_0, metadata, data_dtype), d_0)

  def testAffineQuantizeAllEqual(self):
    d = np.ones(5, dtype=np.float32)
    q, metadata = quantization.quantize_weights(d, np.uint8)
    assert 'scale' in metadata and 'min' in metadata
    self.assertEqual(metadata['scale'], 1.0)
    self.assertEqual(q.dtype, np.uint8)

    de_q = quantization.dequantize_weights(q, metadata, np.float32)
    np.testing.assert_array_equal(de_q, d)

  def testFloatQuantizeAllEqual(self):
    d = np.ones(5, dtype=np.float32)
    q, metadata = quantization.quantize_weights(d, np.float16)
    self.assertDictEqual(metadata, {})

    self.assertEqual(q.dtype, np.float16)
    de_q = quantization.dequantize_weights(q, metadata, np.float32)
    np.testing.assert_array_equal(de_q, d)

  def testAffineQuantizeNegativeFloats(self):
    self._runQuantizeTest(
        -3, -1, np.float32, np.uint8,
        expected_metadata={'scale': 2/255})
    self._runQuantizeTest(
        -3, -1, np.float32, np.uint16,
        expected_metadata={'scale': 2/65536})

  def testAffineQuantizeNegativeAndZeroFloats(self):
    self._runQuantizeTest(
        -3, 0, np.float32, np.uint8,
        expected_metadata={'scale': 3/255})
    self._runQuantizeTest(
        -3, 0, np.float32, np.uint16,
        expected_metadata={'scale': 3/65536})

  def testAffineQuantizeNegativeAndPositiveFloats(self):
    self._runQuantizeTest(
        -3, 3, np.float32, np.uint8,
        expected_metadata={'scale': 6/255})
    self._runQuantizeTest(
        -3, 3, np.float32, np.uint16,
        expected_metadata={'scale': 6/65536})

  def testAffineQuantizeZeroAndPositiveFloats(self):
    self._runQuantizeTest(
        0, 3, np.float32, np.uint8,
        expected_metadata={'scale': 3/255})
    self._runQuantizeTest(
        0, 3, np.float32, np.uint16,
        expected_metadata={'scale': 3/65536})

  def testAffineQuantizePositiveFloats(self):
    self._runQuantizeTest(
        1, 3, np.float32, np.uint8,
        expected_metadata={'scale': 2/255})
    self._runQuantizeTest(
        1, 3, np.float32, np.uint16,
        expected_metadata={'scale': 2/65536})

  def testAffineQuantizeNormalizedFloats(self):
    data = np.array(
        [-0.29098126, -0.24776903, -0.27248842, 0.23848203], dtype=np.float32)
    q, metadata = quantization.quantize_weights(data, np.uint16)
    de_q = quantization.dequantize_weights(q, metadata, data.dtype)
    np.testing.assert_array_almost_equal(de_q, data, decimal=5)

  def testAffineQuantizeNegativeInts(self):
    self._runQuantizeTest(
        -3, -1, np.int32, np.uint8,
        expected_metadata={'scale': 2/255})
    self._runQuantizeTest(
        -3, -1, np.int32, np.uint16,
        expected_metadata={'scale': 2/65536})

  def testAffineQuantizeNegativeAndZeroInts(self):
    self._runQuantizeTest(
        -3, 0, np.int32, np.uint8,
        expected_metadata={'scale': 3/255})
    self._runQuantizeTest(
        -3, 0, np.int32, np.uint16,
        expected_metadata={'scale': 3/65536})

  def testAffineQuantizeNegativeAndPositiveInts(self):
    self._runQuantizeTest(
        -3, 3, np.int32, np.uint8,
        expected_metadata={'scale': 6/255})
    self._runQuantizeTest(
        -3, 3, np.int32, np.uint16,
        expected_metadata={'scale': 6/65536})

  def testAffineQuantizeZeroAndPositiveInts(self):
    self._runQuantizeTest(
        0, 3, np.int32, np.uint8,
        expected_metadata={'scale': 3/255})
    self._runQuantizeTest(
        0, 3, np.int32, np.uint16,
        expected_metadata={'scale': 3/65536})

  def testAffineQuantizePositiveInts(self):
    self._runQuantizeTest(
        1, 3, np.int32, np.uint8,
        expected_metadata={'scale': 2/255})
    self._runQuantizeTest(
        1, 3, np.int32, np.uint16,
        expected_metadata={'scale': 2/65536})

  def testInvalidQuantizationTypes(self):
    # Invalid quantization type
    with self.assertRaises(ValueError):
      quantization.quantize_weights(np.array([]), np.bool)

    # Invalid data dtype for float16 quantization
    with self.assertRaises(ValueError):
      d = np.ones(1, dtype=np.int32)
      quantization.quantize_weights(d, np.float16)

  def testInvalidDequantizationTypes(self):
    # Invalid metadata for affine quantization
    with self.assertRaises(ValueError):
      d = np.ones(1, dtype=np.uint8)
      quantization.dequantize_weights(np.array([]), {})

    # Invalid target dtype for float16 quantization
    with self.assertRaises(ValueError):
      d = np.ones(1, dtype=np.float16)
      quantization.dequantize_weights(d, {}, np.int32)

    # Invalid dequantization type
    with self.assertRaises(ValueError):
      d = np.ones(1, dtype=np.bool)
      quantization.dequantize_weights(d, {})

  def testMapLayerFallthrough(self):
    names = ['conv/0/weight', 'conv/0/bias', 'conv/1/weight', 'conv/1/bias']
    quantization_dtype_map = {'float16': ['conv/0/*'], 'uint8': True}
    dtype_map = quantization.map_layers_to_quantization_dtype(
        names, quantization_dtype_map)

    self.assertDictEqual(dtype_map, {
        'conv/0/weight': np.float16,
        'conv/0/bias': np.float16,
        'conv/1/weight': np.uint8,
        'conv/1/bias': np.uint8
    })

  def testMapLayerConflictingMap(self):
    names = ['conv/0/weight', 'conv/0/bias', 'conv/1/weight', 'conv/1/bias']
    quantization_dtype_map = {'float16': ['conv/0/*'], 'uint8': ['conv/0/bias']}

    with self.assertRaises(ValueError):
      quantization.map_layers_to_quantization_dtype(
          names, quantization_dtype_map)

  def testMapLayerStringToList(self):
    names = ['conv/0/weight', 'conv/0/bias', 'conv/1/weight', 'conv/1/bias']
    quantization_dtype_map = {'float16': '*'}

    dtype_map = quantization.map_layers_to_quantization_dtype(
        names, quantization_dtype_map)

    self.assertDictEqual(dtype_map, {
        'conv/0/weight': np.float16,
        'conv/0/bias': np.float16,
        'conv/1/weight': np.float16,
        'conv/1/bias': np.float16
    })

  def testMapLayerNoDtypeMap(self):
    names = ['conv/0/weight', 'conv/0/bias', 'conv/1/weight', 'conv/1/bias']
    quantization_dtype_map = {}
    dtype_map = quantization.map_layers_to_quantization_dtype(
        names, quantization_dtype_map)

    self.assertDictEqual(dtype_map, {})

  def testMapLayerExactlyOneFallthrough(self):
    names = ['conv/0/weight', 'conv/0/bias', 'conv/1/weight', 'conv/1/bias']
    quantization_dtype_map = {'float16': True, 'uint8': True}

    with self.assertRaises(ValueError):
      quantization.map_layers_to_quantization_dtype(
          names, quantization_dtype_map)

if __name__ == '__main__':
  unittest.main()


text:
import unittest
import os
from pathlib import Path
from citi import Record
import numpy.testing as npt

class TestReadWVIRecord(unittest.TestCase):

    @staticmethod
    def __get_display_memory_filename() -> str:
        relative_path = os.path.join('.', '..', '..', '..')
        this_dir = os.path.dirname(Path(__file__).absolute())
        absolute_path = os.path.join('tests', 'regression_files')
        filename = 'wvi_file.cti'
        return os.path.join(
            this_dir, relative_path, absolute_path, filename
        )

    def setUp(self):
        self.record = Record(self.__get_display_memory_filename())

    def test_file_exists(self):
        os.path.isfile(self.__get_display_memory_filename())

    def test_version(self):
        self.assertEqual(self.record.version, "A.01.01")

    def test_name(self):
        self.assertEqual(self.record.name, "Antonly001")

    def test_comments(self):
        self.assertEqual(len(self.record.comments), 6)
        self.assertEqual(self.record.comments, [
            'SOURCE: 10095059066467',
            'DATE: Fri, Jan 18, 2019, 14:14:44',
            'ANTPOS_TX: 28.4E-3 0E+0 -16E-3 90 270 0',
            'ANTPOS_RX: 28.4E-3 0E+0 -16E-3 90 270 0',
            'ANT_TX: NAH_003',
            'ANT_RX: NAH_003',
        ])

    def test_devices(self):
        self.assertEqual(len(self.record.devices), 0)

    def test_independent_variable(self):
        self.assertEqual(self.record.independent_variable[0], "Freq")
        self.assertEqual(self.record.independent_variable[1], "MAG")
        npt.assert_array_almost_equal(
            self.record.independent_variable[2],
            [100000000.0, 200000000.0]
        )

    def test_data(self):
        self.assertEqual(len(self.record.data), 1)
        self.assertEqual(self.record.data[0][0], 'S11')
        self.assertEqual(self.record.data[0][1], 'RI')
        npt.assert_array_almost_equal(
            self.record.data[0][2],
            [
                complex(8.609423041343689E-1, 4.5087423920631409E-1),
                complex(-6.1961996555328369E-1, -7.2456854581832886E-1),
            ]
        )

max_stars_repo_name: Wave-View-Imaging/citi

max_issues_repo_name: Wave-View-Imaging/citi

max_forks_repo_name: Wave-View-Imaging/citi


text:
import numpy as np
from typing import List, Tuple

from mph import GradedMatrix, groebner_bases, presentation_FIrep

def choose_graded_subbasis(matrix: List[List[Tuple[int, int]]], column_grades: List[List[int]], row_grades: List[List[int]]):
	dense_matrix = np.zeros(shape=(len(row_grades), len(column_grades)), dtype=np.int32)
	for i in range(len(matrix)):
		for entry in matrix[i]:
			dense_matrix[entry[1], i] = entry[0]
	transform = np.random.choice([0, 1], size=(dense_matrix.shape[1], np.random.randint(1, dense_matrix.shape[1]-1) if dense_matrix.shape[1]>2 else 1))
	column_grades_ret = []
	for i in range(transform.shape[1]):
		grade = [0]*len(column_grades[0])
		for j in range(transform.shape[0]):
			if transform[j, i] != 0:
				grade = [max(grade[k], column_grades[j][k]) for k in range(len(grade))]
		column_grades_ret.append(grade)
	return np.matmul(dense_matrix, transform) % 2, column_grades_ret

def generate_rivet_output(high_matrix, column_grades_h, low_matrix, column_grades_l):
	with open("input_file", "w") as f:
		out = "firep\nparameter 1\nparameter 2\n" + str(high_matrix.shape[1]) +" " + str(low_matrix.shape[1]) + " " + str(low_matrix.shape[0])+"\n"
		for column_index in range(high_matrix.shape[1]):
			out += str(column_grades_h[column_index][0]) +" "+str(column_grades_h[column_index][1])+" ; "
			for row_index in range(high_matrix.shape[0]):
				if high_matrix[row_index, column_index] != 0:
					out += str(row_index)+" "
			out += "\n"
		for column_index in range(low_matrix.shape[1]):
			out += str(column_grades_l[column_index][0]) +" "+str(column_grades_l[column_index][1])+" ; "
			for row_index in range(low_matrix.shape[0]):
				if low_matrix[row_index, column_index] != 0:
					out += str(row_index)+" "
			out += "\n"
		f.write(out)

def rivet_benchmark(m: int, n: int, n_parameters: int=2, density: float=0.5, grade_range:int=1000, log_level:str="silent"):
	# Create a random matrix as the low_matrix in the FIrep
	low_matrix = np.random.choice([0, 1], size=(m, n), p=[density, 1-density])
	# Choose random grades for the columns
	column_grades_l = [ np.random.choice(grade_range, size=(n_parameters,)) for _ in range(n) ]
	# Compute the kernel of this graded matrix
	output = groebner_bases(low_matrix, column_grades_l, log_level=log_level)
	if len(output[1].column_grades) == 0:
		return 0
	# Choose a random subbasis of the kernel as the high_matrix in the FIrep
	high_matrix, column_grades_h = choose_graded_subbasis(output[1].matrix, output[1].column_grades, output[1].row_grades)

	generate_rivet_output(high_matrix, column_grades_h, low_matrix, column_grades_l)

	ret = presentation_FIrep(high_matrix, column_grades_h, low_matrix, column_grades_l, log_level=log_level)

def random_FIrep_presentation(m: int, n: int, n_parameters: int=3, density: float=0.5, grade_range:int=100):
	# Create a random matrix as the low_matrix in the FIrep
	low_matrix = np.random.choice([0, 1], size=(m, n), p=[density, 1-density])
	# Choose random grades for the columns
	column_grades_l = [ np.random.choice(grade_range, size=(n_parameters,)) for _ in range(n) ]
	# Compute the kernel of this graded matrix
	output = groebner_bases(low_matrix, column_grades_l)
	if len(output[1].column_grades) == 0:
		return 0
	# Choose a random subbasis of the kernel as the high_matrix in the FIrep
	high_matrix, column_grades_h = choose_graded_subbasis(output[1].matrix, output[1].column_grades, output[1].row_grades)
	return presentation_FIrep(high_matrix, column_grades_h, low_matrix, column_grades_l)

def random_map_gbs(m:int, n:int, n_parameters: int=3, density:float=0.5, grade_range:int=100):
	random_matrix = np.random.choice([0, 1], size=(m, n), p=[density, 1-density])
	# Choose random grades for the columns
	column_grades = [ np.random.choice(grade_range, size=(n_parameters,)) for _ in range(n) ]
	# Compute Groebner bases for the image and kernel of this graded matrix
	image_gb, kernel_gb = groebner_bases(random_matrix, column_grades)
	sparse_random_matrix = [[] for _ in range(random_matrix.shape[1])]
	non_zero_entries = np.nonzero(random_matrix)
	for entry in list(zip(non_zero_entries[1], non_zero_entries[0])):
		sparse_random_matrix[entry[0]].append((random_matrix[entry[1], entry[0]], entry[1]))
	return GradedMatrix(sparse_random_matrix, column_grades, row_grades = [ [ 0 for _ in range(len(column_grades[0])) ] for _ in range(random_matrix.shape[0])]), image_gb, kernel_gb

max_stars_repo_name: MBender/muphasa

max_issues_repo_name: MBender/muphasa

max_forks_repo_name: MBender/muphasa


text:
#!/usr/bin/env python3
import gc
import os
import time
import math
import overpy
import socket
import requests
import threading
import numpy as np
# setup logging
import logging
import logging.handlers
import selfdrive.crash as crash
from common.params import Params
from collections import defaultdict
import cereal.messaging as messaging
from common.realtime import set_realtime_priority

from selfdrive.version import version, dirty

# OPKR 코드 참조.

class MapsdThread(threading.Thread):
    def __init__(self, threadID, name):
        threading.Thread.__init__(self)
        self.threadID =  threadID
        self.name  =  name
        self.logger = logging.getLogger( self.name )

        self.params = Params()
        self.params.put("OpkrMapEnable", "0")
        self.map_enabled = 0
        self.old_map_enable = 0
        self.target_speed_map_counter2 = 0
        self.map_exec_status= False
        self.programRun = True

    # 1,3  map top
    # 2,   backgrand
    # 0,   kill
    def opkr_map_status_read(self):
        self.map_enabled = int(self.params.get("OpkrMapEnable"))

        if self.map_enabled == 2 and self.target_speed_map_counter2 > 0:
            self.target_speed_map_counter2 -= 1
            print( " target_speed_map_counter2 = {}".format( self.target_speed_map_counter2  ))
            if self.target_speed_map_counter2  == 0:
                print( "am start --activity-task-on-home com.opkr.maphack/com.opkr.maphack.MainActivity" )
                os.system("am start --activity-task-on-home com.opkr.maphack/com.opkr.maphack.MainActivity")
                self.programRun = False
                return

        if( self.map_enabled == self.old_map_enable ):
            return

        self.old_map_enable = self.map_enabled

        if self.map_enabled == 0:
            self.map_exec_status = False
            os.system("pkill com.mnsoft.mappyobn")
        elif self.map_exec_status == False:
            self.map_exec_status = True
            os.system("am start com.mnsoft.mappyobn/com.mnsoft.mappy.MainActivity &")  # map 실행.
            if self.map_enabled == 2:
                self.target_speed_map_counter2 = 3   # map 실행후 3초후 backgrand로 전환합니다.
        elif self.map_enabled == 2:  # map backgrand
            os.system("am start --activity-task-on-home com.opkr.maphack/com.opkr.maphack.MainActivity")
        elif self.map_enabled == 3:
            # map return
            os.system("am start --activity-task-on-home com.mnsoft.mappyobn/com.mnsoft.mappy.MainActivity")

    def run(self):
        # OPKR
        self.navi_on_boot = self.params.get_bool("OpkrRunNaviOnBoot")

        if self.navi_on_boot:
            self.params.put("OpkrMapEnable", "2")

        start = time.time()
        while True:
            if time.time() - start < 1.0:
              time.sleep(1.0)
              continue
            else:
              start = time.time()
              self.opkr_map_status_read()

def main():
    #gc.disable()
    #set_realtime_priority(1)
    params = Params()
    dongle_id = params.get("DongleId")

    mt = MapsdThread(1, "/data/media/0/videos/MapsdThread")
    mt.start()

if __name__ == "__main__":
    main()