text:
import itertools
import numpy as np
import cantera as ct
import pandas as pd
import re
import pickle
from .. import simulation as sim
from ...cti_core import cti_processor as ctp

class shockTube(sim.Simulation):

    def __init__(self,pressure:float,temperature:float,observables:list,
                 kineticSens:int,physicalSens:int,conditions:dict,
                 initialTime,finalTime,thermalBoundary,mechanicalBoundary,
                 processor:ctp.Processor=None,cti_path="",save_timeHistories:int=0,
                 save_physSensHistories=0,moleFractionObservables:list=[],
                 absorbanceObservables:list=[],concentrationObservables:list=[],
                 fullParsedYamlFile:dict={},
                 time_shift_value:float = 0, atol:float=1e-15, rtol:float=1e-9,
                 rtol_sens:float=0.0001,
                 atol_sens:float=1e-6):
        '''
        Child class pertaining to shock tube simulations. Inherits all
        attributes and methods from simulations class including __init__().
        Also has its own internal init method due to additional data
        requirements.

        Parameters
        ----------
        pressure : float
            Pressure in [atm].
        temperature : float
            Temperature in [K].
        observables : list
            Species which sensitivity analysis is performed for.
        kineticSens : int
            0 for off, 1 for on.
        physicalSens : int
            0 for off, 1 for on.
        conditions : dict
            Initial mole fractions for species in simulation.
        initialTime : float
            Time to begin simulation from (s).
        finalTime : float
            Time to end simulation (s).
        thermalBoundary : str
            Thermal boundary condition inside the reactor. Shock tubes can
            either be adiabatic or isothermal.
        mechanicalBoundary : TYPE
            Mechanical bondary condition inside the reactor. Shock tubes can
            either be constant pressure or constant volume.
        processor : ctp.Processor, optional
            Loaded cti file. The default is None.
        cti_path : str, optional
            Path of cti file for running. If processor is provided this is not
            needed. The default is "".
        save_timeHistories : int, optional
            Boolean variable describing if time histories for simulation runs
            are saved. 0 for not saved, 1 for saved. The default is 0.
        save_physSensHistories : TYPE, optional
            Boolean variable describing if physical sensitivity time histories
            are saved. 0 for not saved, 1 for saved. The default is 0.
        moleFractionObservables : list, optional
            Species for which experimental data in the form of mole fraction
            time histories will be provided for optimization.
            Kinetic sensitivities are calculated for all these species.
            The default is [].
        absorbanceObservables : list, optional
            Species for which experimental data in the form of summed
            absorption time histories (from some or all of the species) will be
            provided for optimization.
            Kinetic sensitivities are calculated for all these species.
            The default is [].
        concentrationObservables : list, optional
            Species for which experimental data in the form of concentration
            time histories will be provided for optimization.
            Kinetic sensitivities are calculated for all these species.
            The default is [].
        fullParsedYamlFile : dict, optional
            Full dictionary from the parsed shock tube yaml file.
            The default is {}.
        time_shift_value : float, optional
            The numerical value by which the time vector of the simulation
            is shifted in seconds. The default is 0.
        atol : float, optional
            Get the absolute error tolerance. The default is 1e-15.
        rtol : float, optional
            Get the relative error tolerance. The default is 1e-9.
        rtol_sens : float, optional
            Scalar relative error tolerance for sensitivity. The default is 0.0001.
        atol_sens : float, optional
            Scalar absolute error tolerance for sensitivity. The default is 1e-6.

        Returns
        -------
        None.

        '''

        sim.Simulation.__init__(self,pressure,temperature,observables,kineticSens,physicalSens,
                                conditions,processor,cti_path)
        self.initialTime = initialTime
        self.finalTime = finalTime
        self.thermalBoundary = thermalBoundary
        self.mechanicalBoundary = mechanicalBoundary
        self.kineticSensitivities= None
        self.timeHistory = None
        self.experimentalData = None
        self.concentrationObservables = concentrationObservables
        self.moleFractionObservables = moleFractionObservables
        self.absorbanceObservables = absorbanceObservables
        self.fullParsedYamlFile =  fullParsedYamlFile
        self.time_shift_value = time_shift_value

        if save_timeHistories == 1:
            self.timeHistories=[]
            self.timeHistoryInterpToExperiment = None
            self.pressureAndTemperatureToExperiment = None
        else:
            self.timeHistories=None
        if save_physSensHistories == 1:
            self.physSensHistories = []
        self.setTPX()
        self.dk = [0]
        self.atol=atol
        self.rtol=rtol
        self.rtol_sensitivity=rtol_sens
        self.atol_sensitivity=atol_sens

    def printVars(self):
        '''
        Prints variables associated with the reactor initial conditions.

        '''

        print('initial time: {0}\nfinal time: {1}\n'.format(self.initialTime,self.finalTime),
              '\nthermalBoundary: {0}\nmechanicalBoundary: {1}'.format(self.thermalBoundary,self.mechanicalBoundary),
              '\npressure: {0}\ntemperature: {1}\nobservables: {2}'.format(self.pressure,self.temperature,self.observables),
              '\nkineticSens: {0}\nphysicalSens: {1}'.format(self.kineticSens,self.physicalSens),
              '\nTPX: {0}'.format(self.processor.solution.TPX)
              )
    #maybe unify paths with cti file?, also really fix the python styling
    def write_time_histories(self, path=''):
        if self.timeHistories == None:
            print("Error: this simulation is not saving time histories, reinitialize with flag")
            return -1
        if path=='':
            path = './time_histories.time'
        pickle.dump(self.timeHistories,open(path,'wb'))
        return 0

    def write_physSensHistories(self, path=''):
        if self.physSensHistories == None:
            print("Error: this simulation is not saving time histories, reinitialize with flag")
            return -1
        if path=='':
            path = './physSensHistories.sens'
            pickle.dump(self.physSensHistories,open(path,'wb'))
        return 0

    def load_physSensHistories(self, path=''):
        if self.physSensHistories == None:
            print("Error: this simulation is not saving time histories, reinitialize with flag")
            return -1
        if path=='':
            path = './physSensHistories.sens'
        pickle.load(self.physSensHistories,open(path,'wb'))
        return 0

    #note this is destructive, the original timeHistories are overwritten, run before any runs
    #same write is destructive by default
    def load_time_histories(self, path=''):
        if self.timeHistories == None:
            print("Error: this simulation is not saving time histories, reinitialize with flag")
            return -1
        if path=='':
            path = './time_histories.time'
        pickle.load(self.timeHistories,open(path,'wb'))
        return 0
    def settingShockTubeConditions(self):
        '''
        Determine the mechanical and thermal boundary conditions for a
        shock tube based on what was initialized.
        '''

        #assigning the thermal boundary variable
        if re.match('[aA]diabatic',self.thermalBoundary):
            energy = 'on'
        elif re.match('[iI]sothermal',self.thermalBoundary):
            energy = 'off'
        else:
            raise Exception('Please specify a thermal boundary condition, adiabatic or isothermal')
        #assigning the mehcanical boundary variable
        if re.match('[Cc]onstant [Pp]ressure',self.mechanicalBoundary):
            mechBoundary = 'constant pressure'
        elif re.match('[Cc]onstant [Vv]olume',self.mechanicalBoundary):
            mechBoundary = 'constant volume'
        else:
            raise Exception('Please specifiy a mehcanical boundary condition, constant pressure or constant volume')
        #return the thermal and mechanical boundary of the shock tube
        return energy,mechBoundary

    def sensitivity_adjustment(self,temp_del:float=0.0,
                               pres_del:float=0.0,
                               spec_pair:(str,float)=('',0.0)):
        '''
        Appends the sensitivity adjustment to list, and calls sensitivity adjustment
        function from the simulations class, to adjust P,T,X for the sensitivity
        calculation

        Parameters
        ----------
        temp_del : float, optional
            The decimal value of the percent by which temperature is perturbed.
            The default is 0.0.
        pres_del : float, optional
            The decimal value of the percent by which pressure is perturbed.
            The default is 0.0.
        spec_pair : (str,float), optional
             The string of a species and the decimal value of the percent by
             which that species is perturbed . The default is ('',0.0).

        Returns
        -------
        data : Pandas Data Frame
            Time history of the perturbed simulation.

        '''

        if temp_del != 0.0:
            self.dk.append(temp_del)
        if pres_del != 0.0:
            self.dk.append(pres_del)
        if spec_pair[1] != 0.0:
            self.dk.append(spec_pair[1])

        kin_temp = self.kineticSens
        self.kineticSens = 0
        data = sim.Simulation.sensitivity_adjustment(self,temp_del,pres_del,spec_pair)
        self.kineticSens = kin_temp
        return data

    def run(self,initialTime:float=-1.0, finalTime:float=-1.0):
        '''
        Run the shock tube simulation

        Parameters
        ----------
        initialTime : float, optional
            The time at which the reactor simulation begins, in seconds. The default is -1.0.
        finalTime : float, optional
            The time at which the reactor simulation ends, in seconds. The default is -1.0.

        Returns
        -------
        timeHistory: Pandas DataFrame
            Time history of simulation containing temperature, pressurea and
            species results.
        kineticSensitivities: numpy array
            three dimensional numpy array: (time x reaction x observable).

        '''

        if initialTime == -1.0:
            initialTime = self.initialTime
        if finalTime == -1.0:
            finalTime = self.finalTime
        self.timeHistory = None
        self.kineticSensitivities= None #3D numpy array, columns are reactions with timehistories, depth gives the observable for those histories
        conditions = self.settingShockTubeConditions()
        mechanicalBoundary = conditions[1]
        #same solution for both cp and cv sims
        if mechanicalBoundary == 'constant pressure':
            shockTube = ct.IdealGasConstPressureReactor(self.processor.solution,
                                                        name = 'R1',
                                                        energy = conditions[0])
        else:
            shockTube = ct.IdealGasReactor(self.processor.solution,
                                           name = 'R1',
                                           energy = conditions[0])
        sim = ct.ReactorNet([shockTube])
        sim.rtol=self.rtol
        sim.atol=self.atol
        #print(sim.rtol_sensitivity,sim.atol_sensitivity)
        sim.rtol_sensitivity=self.rtol_sensitivity
        sim.atol_sensitivity=self.atol_sensitivity

        columnNames = [shockTube.component_name(item) for item in range(shockTube.n_vars)]
        columnNames = ['time']+['pressure']+columnNames
        self.timeHistory = pd.DataFrame(columns=columnNames)

        if self.kineticSens == 1:
            for i in range(self.processor.solution.n_reactions):
                shockTube.add_sensitivity_reaction(i)
            dfs = [pd.DataFrame() for x in range(len(self.observables))]
            tempArray = [np.zeros(self.processor.solution.n_reactions) for x in range(len(self.observables))]

        t = self.initialTime
        counter = 0
        #print(sim.rtol_sensitivity,sim.atol_sensitivity)
        while t < self.finalTime:
            t = sim.step()
            if mechanicalBoundary =='constant volume':
                state = np.hstack([t,shockTube.thermo.P,shockTube.mass,shockTube.volume,
                               shockTube.T, shockTube.thermo.X])
            else:
                state = np.hstack([t,shockTube.thermo.P, shockTube.mass,
                               shockTube.T, shockTube.thermo.X])

            self.timeHistory.loc[counter] = state
            if self.kineticSens == 1:
                counter_1 = 0
                for observable,reaction in itertools.product(self.observables, range(self.processor.solution.n_reactions)):
                    tempArray[self.observables.index(observable)][reaction] = sim.sensitivity(observable,
                                                                                reaction)
                    counter_1 +=1
                    if counter_1 % self.processor.solution.n_reactions == 0:
                        dfs[self.observables.index(observable)] = dfs[self.observables.index(observable)].append(((
                            pd.DataFrame(tempArray[self.observables.index(observable)])).transpose()),
                            ignore_index=True)
            counter+=1

        if self.timeHistories != None:

            self.timeHistory.time = self.timeHistory.time + self.time_shift_value

            #self.timeHistory.time = self.timeHistory.time + 0

            self.timeHistories.append(self.timeHistory)
            ############################################################

        if self.kineticSens == 1:
            numpyMatrixsksens = [dfs[dataframe].values for dataframe in range(len(dfs))]
            self.kineticSensitivities = np.dstack(numpyMatrixsksens)
            return self.timeHistory,self.kineticSensitivities
        else:
            return self.timeHistory

    #interpolate the most recent time history against the oldest by default
    #working_data used if have list not pandas frame
    #return more data about what was interpolated in a tuple?
    def interpolate_time(self,index:int=None,time_history=None):
        '''
        This function interpolates and returns the most recent time history
        against the original time history by default, unless a specific time
        history index  or time history is passed in. If an index is passed in
        then an interpolated time history associated with that index in the
        list is returned. If a specific time_history is passed in then an
        interpolated version of that time history is returned.

        Parameters
        ----------
        index : int, optional
            The index value of the specific time history to be interpolated.
            The default is None.
        time_history : Pandas Data Frame, optional
            Pandas dataframe containing the time history from a simulation.
            The default is None.

        Returns
        -------
        interpolatedTimeHistory: Pandas Data Frame
            Interpolated time history.

        '''

        if self.timeHistories == None:
            print("Error: this simulation is not saving time histories, reinitialize with flag")
            return -1
        else:
            if index is not None and time_history is not None:
                print("Error: can only specify one of index, time_history")
                return -1
            if index is None:
                index = -1

            if time_history is None:
                return self.interpolation(self.timeHistories[0],self.timeHistories[index],["temperature","pressure"]+self.observables)
            else:
                return self.interpolation(self.timeHistories[0],time_history,["temperature","pressure"]+self.observables)

    #assumes most recent time histories are the correct ones to interpolate on
    #interpolates agains the original time history
    def interpolate_species_adjustment(self):
        '''
        This function interpolates the time history of a species adjustment run,
        against the original time history.
        '''
        interpolated_data = []
        species_to_loop = set(self.conditions.keys()).difference(['Ar','AR','HE','He','Kr','KR','Xe','XE','NE','Ne'])

        for x in range(0,len(species_to_loop)):
            interpolated_data.insert(0,self.interpolate_time(index=-1-x))

        return interpolated_data
    def interpolate_species_sensitivities(self):
        '''
        This function interpolates the time history of a species sensitivity
        against the original time history.
        '''
        interpolated_data = self.interpolate_species_adjustment()
        interpolated_sens = []

        ind_off = len(self.timeHistories)-len(set(self.conditions.keys()).difference(['Ar','AR','HE','He','Kr','KR','Xe','XE','NE','Ne']))
        for i,th in enumerate(interpolated_data):
            ind = ind_off + i

            interpolated_sens.append(self.interpolate_physical_sensitivities(index=ind,time_history=th))

        return interpolated_sens

    #interpolate a range of time histories agains the original
    #possibly add experimental flag to do range with exp data
    #end is exclusive
    #start here tomorrow
    def interpolate_range(self,begin:int,end:int):
        '''
        Function that defines a time range for interpolation.

        Parameters
        ----------
        begin : int
            Time to begin interpolation.
        end : int
            Time to end interpolation.

        Returns
        -------
        interpolated_data : Pandas Data Frame
            Interpolated time history.

        '''
        if begin<0 or end>len(self.timeHistories):
            print("Error: invalid indices")
        if self.timeHistories == None:
            print("Error: simulation is not saving time histories")

        interpolated_data = []
        for x in range(begin,end):
            interpolated_data.append(self.interpolate_time(index=x))
        return interpolated_data

    #interpolates agains the original time history
    def interpolate_physical_sensitivities(self, index:int=-1,
                                           time_history=None):
        '''
        This function interpolates the time history of a physical sensitivity,
        excluding species, against the original time history and returns a
        numpy array.

        Parameters
        ----------
        index : int, optional
            The index value of the specific time history to be interpolated.
            The default is -1.
        time_history : Pandas Data Frame, optional
            Pandas dataframe containing the time history from a simulation.
            The default is None.

        Returns
        -------
        sensitivity: Pandas Data Frame.
        Physical sensitivity Data Frame.

        '''

        interpolated_time = self.interpolate_time(index) if time_history is None else time_history
        #print(interpolated_time)
        #calculate which dk
        dk = self.dk[index]
        sensitivity = self.sensitivityCalculation(self.timeHistories[0][self.observables],
                                                  interpolated_time[self.observables],self.observables,dk)
        if self.physSensHistories != None:
            self.physSensHistories.append(sensitivity)
#        print('this is sensitivity')
#        print(sensitivity)
        return sensitivity

    #returns a 3D array of interpolated time histories corrosponding to physical sensitivities
    def interpolate_experimental_kinetic(self, pre_interpolated = []):
        '''
        This function interpolates kinetic sensitivities to experimental data.

        Parameters
        ----------
        pre_interpolated : list, optional
            List of kinetic sensitivties to be interpolated. If not provided
            the function will look for kinetic sensitivities that were saved
            as an attribute.
            The default is [].

        Returns
        -------
        flipped: numpy array
            Interpolated 3d array containing sensitivities.

        '''

        if self.experimentalData == None:
            print("Error: experimental data must be loaded")
            return -1
        if len(pre_interpolated) == 0 and not self.kineticSensitivities.any():
            print("Error: must specify pre_interpolated or have kineticSensitivities run first")
            return -1
        elif len(pre_interpolated)>0:
            array = pre_interpolated
        else:
            array = self.kineticSensitivities
        exp_interp_array = []
        #if len(self.experimentalData) < array.shape[2]:
        #    print("Error: mismatch between kineticSensitivities observables and given experimental data")
        #    return -1
        #exp data and kineticSensitivities must match in size and order of observables
        for i,frame in enumerate(self.experimentalData):
            if i > array.shape[2]:
                break
            sheet = array[:,:,i]
            exp_interp_array.append([])
            for time_history in sheet.T:
                # new_history = np.interp(frame.ix[:,0],
                #                        self.timeHistories[0]['time'],
                #                        time_history)
                new_history = np.interp(frame.iloc[:,0],
                             self.timeHistories[0]['time'],
                             time_history)
                new_history = new_history.reshape((new_history.shape[0],
                                                  1))
                exp_interp_array[-1].append(new_history)
        flipped = []
        for x in exp_interp_array:
            flipped.append(np.hstack(x))
        return flipped

    def map_and_interp_ksens(self,time_history=None):
        '''
        This function maps kinetic sensitivity calculations returned from cantera
        to kineitcs parameters A,n and Ea, as well as interpolates them to the
        corresponding experimental data. It returns a dictonary containing the
        interpolated kinetic senstivities.

        Parameters
        ----------
        time_history : Pandas Data Frame, optional
            The original time history is required for obtaining temperature
            values in order to do the mapping. The default is None.

        Returns
        -------
        dict
            Dictionary containing the mapping for kinetic sensitivities.

        '''

        A = self.kineticSensitivities
        N = np.zeros(A.shape)
        Ea = np.zeros(A.shape)
        for i in range(0,A.shape[2]):
            sheetA = A[:,:,i] #sheet for specific observable
            for x,column in enumerate(sheetA.T):
                N[:,x,i]= np.multiply(column,np.log(self.timeHistories[0]['temperature'])) if time_history is None else np.multiply(column,np.log(time_history['temperature']))
                #not sure if this mapping is correct, check with burke and also update absorption mapping
                #to_mult_ea = np.divide(-1,np.multiply(1/ct.gas_constant,self.timeHistories[0]['temperature'])) if time_history is None else np.divide(-1,np.multiply(ct.gas_constant,time_history['temperature']))
                to_mult_ea = np.divide(-1,np.multiply(1,self.timeHistories[0]['temperature'])) if time_history is None else np.divide(-1,np.multiply(1,time_history['temperature']))
                Ea[:,x,i]= np.multiply(column,to_mult_ea)

        return {'A':self.interpolate_experimental_kinetic(A),
                'N':self.interpolate_experimental_kinetic(N),
                'Ea':self.interpolate_experimental_kinetic(Ea)}

    #assumes pre_interpolated has been interpolated against the original time history
    #assumes pre_interpolated is a list of dataframes where each dataframe is a time history
    #single is a single dataframe representing one time history/run of the simulation
    def interpolate_experimental(self,pre_interpolated = [], single = None):
        '''
        This function interpolates the time history of a physical sensitivity
        against the corresponding experimental data and returns a pandas
        dataframe.

        Parameters
        ----------
        pre_interpolated : list, optional
            has been interpolated against the original time history,
        assumes pre_interpolated is a list of dataframes where each dataframe
        is a time history. The default is [].
        single : Pandas Data Frame, optional
            Single Pandas Data Frame to be interpolated against experimental
            data. The default is None.

        Returns
        -------
        int_exp: Pandas Data Frame
            Interpolated data frame.

        '''

        if self.timeHistories == None:
            print("Error: can't interpolate without time histories")
            return -1
        if self.experimentalData == None:
            print("Error: must have experimental data before interpolation")
            return -1
        if len(pre_interpolated)!=0 and single != None:
            print("Error: can only specify one of pre_interpolated, single")
        if single is not  None:
            pre_interpolated = [single]

        int_exp = []
        #This portion of the function removes the pressure and the temperature from the pre_interpolated frames
        #so that when it gets interpolated against experimental data the temperature and pressure are not there
        if isinstance(pre_interpolated[0],pd.DataFrame):
            if 'pressure' in pre_interpolated[0].columns.tolist() and 'temperature' in pre_interpolated[0].columns.tolist():
                pre_interpolated = [df.drop(columns=['temperature','pressure']) for df in pre_interpolated]

        #make sure you put the observables list in the correct order
        #check what order the experimental list is parsed in
        #making new observables list
        if single is not None:
            mole_fraction_and_concentration_observables= self.moleFractionObservables + self.concentrationObservables

            mole_fraction_and_concentration_observables   = [x for x in mole_fraction_and_concentration_observables if x is not None]
            for time_history in pre_interpolated:
                array_list = []
                max_size = 0
                for i, observable in enumerate(mole_fraction_and_concentration_observables):
                    interpolated_column = np.interp(self.experimentalData[i]['Time'].values,
                                                    self.timeHistories[0]['time'].values,
                                                    time_history[observable].values)

                    interpolated_column = np.reshape(interpolated_column,((interpolated_column.shape[0],1)))
                    array_list.append(interpolated_column)
                    max_size = max(interpolated_column.shape[0],max_size)
                padded_arrays = []
                for arr in array_list:
                    if arr.shape[0] < max_size:
                        padded_arrays.append(np.pad(arr,
                                                    ((0,max_size-arr.shape[0]),(0,0)),
                                                    'constant',constant_values = np.nan))
                    else:
                        padded_arrays.append(arr)

                np_array = np.hstack((padded_arrays))
                new_frame = pd.DataFrame(np_array)
                int_exp.append(new_frame)

            for x in int_exp:
                x.columns = mole_fraction_and_concentration_observables

            return int_exp[0]

        #check and make sure this part actually works for what we want  for interpolating the pre interpolated time histories
        #make sure we are getting the correct columns
        else:
            for time_history in pre_interpolated:
                array_list = []
                max_size = 0
                for i,frame in enumerate(self.experimentalData): #each frame is data for one observable
                    if i>len(self.observables):
                        break
                    #change these bboth to use observable
                    # interpolated_column= np.interp(frame.ix[:,0],
                    #                                self.timeHistories[0]['time'],
                    #                                time_history.ix[:,i])

                    interpolated_column= np.interp(frame.iloc[:,0],
                                                   self.timeHistories[0]['time'],
                                                   time_history.iloc[:,i])

                    interpolated_column= np.reshape(interpolated_column,
                                                ((interpolated_column.shape[0],1)))
                    array_list.append(interpolated_column)
                    max_size = max(interpolated_column.shape[0],max_size)
                padded_arrays= []
                for arr in array_list:
                    if arr.shape[0] < max_size:
                        padded_arrays.append(np.pad(arr,
                                            ((0,max_size - arr.shape[0]),(0,0)),
                                            'constant',constant_values=np.nan))
                    else:
                        padded_arrays.append(arr)
                np_array = np.hstack((padded_arrays))
                new_frame = pd.DataFrame(np_array)
                int_exp.append(new_frame)

            for x in int_exp:
                x.columns = self.observables[0:len(self.experimentalData)]
            if single is not  None:
                return int_exp[0]
            else:
                return int_exp

    def interpolation(self,originalValues,newValues, thingBeingInterpolated):
        '''
        This function is the base interpolation function, interpolating one set
        of data to another on a per time basis and returns a pandas dataframe.

        Parameters
        ----------
        originalValues : Pandas Data Frame
            Original dataframe of time history.
        newValues : Pandas Data Frame
            New dataframe of time history.
        thingBeingInterpolated : list
            List of observable names.

        Returns
        -------
        interpolatedData: Pandas Data Frame
            Interpolated Data Frame.

        '''
        #interpolating time histories to original time history

        if isinstance(originalValues,pd.DataFrame) and isinstance(newValues,pd.DataFrame):
            tempDfForInterpolation = newValues[thingBeingInterpolated]
            #tempListForInterpolation = [tempDfForInterpolation.ix[:,x].values for x in range(tempDfForInterpolation.shape[1])]
            tempListForInterpolation = [tempDfForInterpolation.iloc[:,x].values for x in range(tempDfForInterpolation.shape[1])]
            interpolatedData = [np.interp(originalValues['time'].values,newValues['time'].values,tempListForInterpolation[x]) for x in range(len(tempListForInterpolation))]
            interpolatedData = [pd.DataFrame(interpolatedData[x]) for x in range(len(interpolatedData))]
            interpolatedData = pd.concat(interpolatedData, axis=1,ignore_index=True)
            interpolatedData.columns = thingBeingInterpolated
        else:
            print("Error: values must be pandas dataframes")
            return -1
        return interpolatedData

    def sensitivityCalculation(self,originalValues,newValues,thingToFindSensitivtyOf,dk=.01):
        '''
         This function calculates log/log sensitivity and returns a
         pandas dataframe.

        Parameters
        ----------
        originalValues : Pandas Data Frame
            Original dataframe of time history.
        newValues : Pandas Data Frame
             New dataframe of time history.
        thingToFindSensitivtyOf : list
            List of observable names.
        dk : float, optional
            The decimal value of the percentage by which the original value
            was perturbed. The default is .01.

        Returns
        -------
        sensitivity: Pandas Data Frame
            Sensitivity of observables.

        '''

        if isinstance(originalValues,pd.DataFrame) and isinstance(newValues,pd.DataFrame):
            newValues.columns = thingToFindSensitivtyOf
            newValues = newValues.applymap(np.log)
            originalValues = originalValues.applymap(np.log)
            #tab

            sensitivity = (newValues.subtract(originalValues)/dk)
            return sensitivity
        else:
            print("Error: wrong datatype, both must be pandas data frames")
            return -1

    def importExperimentalData(self,csvFileList):
        '''
        This function imports experimental data in csv format and returns a
        list of pandas dataframes.

        Parameters
        ----------
        csvFileList : list
            List of csv file directories.

        Returns
        -------
        experimentalData : list
            Experimental data from csv files as pandas data frames
            stored in a list.

        '''

        print('Importing shock tube data the following csv files...')
        print(csvFileList)
        experimentalData = [pd.read_csv(csv) for csv in csvFileList]
        experimentalData = [experimentalData[x].dropna(how='any') for x in range(len(experimentalData))]
        experimentalData = [experimentalData[x].apply(pd.to_numeric, errors = 'coerce').dropna() for x in range(len(experimentalData))]
        for x in range(len(experimentalData)):
            experimentalData[x] = experimentalData[x][~(experimentalData[x][experimentalData[x].columns[1]] < 0)]
        self.experimentalData = experimentalData
        return experimentalData

    def savingInterpTimeHistoryAgainstExp(self,timeHistory):
        '''
        This function writes the time history which is interpolated against
        experimetnal data as a class object.

        Parameters
        ----------
        timeHistory : Pandas Data Frame
            Time history interpolated against experimental data .

        Returns
        -------
        timeHistoryInterpToExperiment: Pandas Data Frame
            Time history that is being saved.

        '''

        self.timeHistoryInterpToExperiment = timeHistory

    def interpolatePressureandTempToExperiment(self,simulation,experimental_data):
        '''
        This function interpolates the pressure and temperature time history
        from a simulation against the corresponding experimental data.

        Parameters
        ----------
        simulation : class variable
             The simulation assoicated with the time history to be
             interpolated.
        experimental_data : list
            List of corresponding experimental data dataframtes .

        Returns
        -------
        list_of_df : list
            List of interpolated data frames for pressure and temperature.

        '''

        p_and_t = ['pressure','temperature']
        list_of_df = []
        for df in experimental_data:
            temp = []
            for variable in p_and_t:
                interpolated_data = np.interp(df['Time'],simulation.timeHistories[0]['time'],simulation.timeHistories[0][variable])
                interpolated_data = interpolated_data.reshape((interpolated_data.shape[0],1))
                temp.append(interpolated_data)
            temp = np.hstack(temp)
            temp = pd.DataFrame(temp)
            temp.columns = p_and_t
            list_of_df.append(temp)
            self.pressureAndTemperatureToExperiment = list_of_df

        return list_of_df

    def calculate_time_shift_sensitivity(self,simulation,experimental_data,dk=1e-8):
        '''
        This function interpolates the pressure and temperature time history
        from a simulation against the corresponding experimental data.

        Parameters
        ----------
        simulation : class variable
            The simulation assoicated with the time history to
            be interpolated.
        experimental_data : list
             List of corresponding experimental data dataframtes .
        dk : float, optional
            Decimal percentage by which time values are perturbed. Default is
            1e-8.

        Returns
        -------
        time_shift_sensitivity : Pandas Data Frame
            List of Pandas Data Frames containing time shift sensitivity
            values.

        '''

        lst_obs = simulation.moleFractionObservables + simulation.concentrationObservables
        lst_obs = [i for i in lst_obs if i]

        one_percent_of_average = dk

        original_time = simulation.timeHistories[0]['time']
        new_time = original_time + one_percent_of_average

        #interpolate to the orignal time
        interpolated_against_original_time = []
        for i,obs in enumerate(lst_obs):
            interpolated_original_observable_against_original_time = np.interp(original_time,new_time,simulation.timeHistories[0][lst_obs[i]])
            s1 = pd.Series(interpolated_original_observable_against_original_time,name=lst_obs[i])
            interpolated_against_original_time.append(s1)

        observables_interpolated_against_original_time_df = pd.concat(interpolated_against_original_time,axis=1)

        #calculate sensitivity

        calculated_sensitivity = []
        for i,obs in enumerate(lst_obs):

           sens = (observables_interpolated_against_original_time_df[obs].apply(np.log) - simulation.timeHistories[0][obs].apply(np.log))/one_percent_of_average
           s1 = pd.Series(sens,name=lst_obs[i])
           calculated_sensitivity.append(s1)

        calculated_sensitivity_df = pd.concat(calculated_sensitivity,axis=1)

        interpolated_sensitivity = []
        for i,obs in enumerate(lst_obs):
            interpolated_sensitivty_per_original_observable = np.interp(experimental_data[i]['Time'],simulation.timeHistories[0]['time'],calculated_sensitivity_df[obs])
            s1 = pd.Series(interpolated_sensitivty_per_original_observable,name=obs)
            interpolated_sensitivity.append(s1)

        interpolated_sensitivity_df = pd.concat(interpolated_sensitivity,axis=1)

        time_shift_sensitivity = interpolated_sensitivity_df

        self.time_shift_sensitivity = time_shift_sensitivity
        average_time=1
        self.average_time = average_time
        return time_shift_sensitivity

max_stars_repo_name: carlylagrotta/MSI

max_issues_repo_name: TheBurkeLab/MSI

max_forks_repo_name: TheBurkeLab/MSI


text:
"""
(c) DI Dominik Hirner BSc.
Institute for graphics and vision (ICG)
University of Technology Graz, Austria
e-mail: dominik.hirner@tugraz.at
"""
import sys
import numpy as np
import cv2
import re
import numpy.matlib
import torch
import torch.nn as nn
from torch.autograd import Variable
from PIL import Image
from typing import Tuple
import torch.nn.functional as F
from guided_filter_pytorch.guided_filter import GuidedFilter
import argparse

def main():

    parser = argparse.ArgumentParser(description='FC-DCNN disparity networ')

    parser.add_argument('--left', help='path to left rectified image')
    parser.add_argument('--weights', help='path to the trained weights')
    parser.add_argument('--right', help='path to right image')
    parser.add_argument('--max_disp', help='disparity search range', type=int)
    parser.add_argument('--out', help='path and/or name for output')
    args = parser.parse_args()

    weight = args.weights
    left_im = args.left
    right_im = args.right
    max_disp = args.max_disp
    out_fn = args.out

    branch.load_state_dict(torch.load(weight))
    disp_s, disp = TestImage(left_im, right_im, max_disp, out_fn, filtered = True, lr_check = True, fill = True)

num_conv_feature_maps = 64
class SiameseBranch(nn.Module):
    def __init__(self,img_ch=3):
        super(SiameseBranch,self).__init__()

        self.Tanh = nn.Tanh()
        self.Conv1 = nn.Conv2d(img_ch, num_conv_feature_maps, kernel_size = 3,stride=1,padding = 1,dilation = 1, bias=True)
        self.Conv2 = nn.Conv2d(num_conv_feature_maps, num_conv_feature_maps, kernel_size=3,stride=1,padding = 1,dilation = 1, bias=True)
        self.Conv3 = nn.Conv2d(2*num_conv_feature_maps, num_conv_feature_maps, kernel_size=3,stride=1,padding = 1,dilation = 1, bias=True)
        self.Conv4 = nn.Conv2d(3*num_conv_feature_maps, num_conv_feature_maps, kernel_size=3,stride=1,padding = 1,dilation = 1,bias=True)
        self.Conv5 = nn.Conv2d(4*num_conv_feature_maps, num_conv_feature_maps, kernel_size=3,stride=1,padding = 1,dilation = 1, bias=True)

    def forward(self,x_in):

        x1 = self.Conv1(x_in)
        x1 = self.Tanh(x1)

        x2 = self.Conv2(x1)
        x2 = self.Tanh(x2)

        d2 = torch.cat((x1,x2),dim=1)

        x3 = self.Conv3(d2)
        x3 = self.Tanh(x3)

        d3 = torch.cat((x1,x2,x3),dim=1)

        x4 = self.Conv4(d3)
        x4 = self.Tanh(x4)

        d4 = torch.cat((x1,x2,x3,x4),dim=1)

        x5 = self.Conv5(d4)

        return x5

Tensor = torch.cuda.FloatTensor
cos = torch.nn.CosineSimilarity()

branch = SiameseBranch()
branch = branch.cuda()

##python3 version!!!!
def readPFM(file):
    file = open(file, 'rb')

    color = None
    width = None
    height = None
    scale = None
    endian = None

    header = file.readline().decode('utf-8').rstrip()
    if header == 'PF':
        color = True
    elif header == 'Pf':
        color = False
    else:
        raise Exception('Not a PFM file.')

    dim_match = re.match(r'^(\d+)\s(\d+)\s$', file.readline().decode('utf-8'))
    if dim_match:
        width, height = map(int, dim_match.groups())
    else:
        raise Exception('Malformed PFM header.')

    scale = float(file.readline().decode('utf-8').rstrip())
    if scale < 0:  # little-endian
        endian = '<'
        scale = -scale
    else:
        endian = '>'  # big-endian

    data = np.fromfile(file, endian + 'f')
    shape = (height, width, 3) if color else (height, width)

    data = np.reshape(data, shape)
    data = np.flipud(data)
    return data, scale

def writePFM(file, image, scale=1):
    file = open(file, 'wb')

    color = None

    if image.dtype.name != 'float32':
        raise Exception('Image dtype must be float32.')

    image = np.flipud(image)

    if len(image.shape) == 3 and image.shape[2] == 3:  # color image
        color = True
    elif len(image.shape) == 2 or len(image.shape) == 3 and image.shape[2] == 1:  # greyscale
        color = False
    else:
        raise Exception('Image must have H x W x 3, H x W x 1 or H x W dimensions.')

    file.write('PF\n'.encode() if color else 'Pf\n'.encode())
    file.write('%d %d\n'.encode() % (image.shape[1], image.shape[0]))

    endian = image.dtype.byteorder

    if endian == '<' or endian == '=' and sys.byteorder == 'little':
        scale = -scale

    file.write('%f\n'.encode() % scale)

    image.tofile(file)

def _compute_binary_kernel(window_size: Tuple[int, int]) -> torch.Tensor:
    r"""Creates a binary kernel to extract the patches. If the window size
    is HxW will create a (H*W)xHxW kernel.
    """
    window_range: int = window_size[0] * window_size[1]
    kernel: torch.Tensor = torch.zeros(window_range, window_range)
    for i in range(window_range):
        kernel[i, i] += 1.0
    return kernel.view(window_range, 1, window_size[0], window_size[1])

def _compute_zero_padding(kernel_size: Tuple[int, int]) -> Tuple[int, int]:
    r"""Utility function that computes zero padding tuple."""
    computed: Tuple[int, ...] = tuple([(k - 1) // 2 for k in kernel_size])
    return computed[0], computed[1]

class MedianBlur(nn.Module):
    r"""Blurs an image using the median filter.

    Args:
        kernel_size (Tuple[int, int]): the blurring kernel size.

    Returns:
        torch.Tensor: the blurred input tensor.

    Shape:
        - Input: :math:`(B, C, H, W)`
        - Output: :math:`(B, C, H, W)`

    Example:
        >>> input = torch.rand(2, 4, 5, 7)
        >>> blur = kornia.filters.MedianBlur((3, 3))
        >>> output = blur(input)  # 2x4x5x7
    """

    def __init__(self, kernel_size: Tuple[int, int]) -> None:
        super(MedianBlur, self).__init__()
        self.kernel: torch.Tensor = _compute_binary_kernel(kernel_size)
        self.padding: Tuple[int, int] = _compute_zero_padding(kernel_size)

    def forward(self, input: torch.Tensor):  # type: ignore
        if not torch.is_tensor(input):
            raise TypeError("Input type is not a torch.Tensor. Got {}"
                            .format(type(input)))
        if not len(input.shape) == 4:
            raise ValueError("Invalid input shape, we expect BxCxHxW. Got: {}"
                             .format(input.shape))
        # prepare kernel
        b, c, h, w = input.shape
        tmp_kernel: torch.Tensor = self.kernel.to(input.device).to(input.dtype)
        kernel: torch.Tensor = tmp_kernel.repeat(c, 1, 1, 1)

        # map the local window to single vector
        features: torch.Tensor = F.conv2d(
            input, kernel, padding=self.padding, stride=1, groups=c)
        features = features.view(b, c, -1, h, w)  # BxCx(K_h * K_w)xHxW

        # compute the median along the feature axis
        median: torch.Tensor = torch.median(features, dim=2)[0]
        return median

# functiona api
def median_blur(input: torch.Tensor,
                kernel_size: Tuple[int, int]) -> torch.Tensor:
    r"""Blurs an image using the median filter.

    See :class:`~kornia.filters.MedianBlur` for details.
    """
    return MedianBlur(kernel_size)(input)

def filterCostVolMedian(cost_vol):

    d,h,w = cost_vol.shape
    cost_vol = cost_vol.unsqueeze(0)

    for disp in range(d):
        cost_vol[:,disp,:,:] = median_blur(cost_vol[:,disp,:,:].unsqueeze(0), (5,5))

    return torch.squeeze(cost_vol)

def filterCostVolBilat(cost_vol,left):

    left = np.mean(left,axis=2)
    leftT = Variable(Tensor(left))
    leftT = leftT.unsqueeze(0).unsqueeze(0)

    d,h,w = cost_vol.shape

    f = GuidedFilter(8,10).cuda() #0.001

    for disp in range(d):
        cur_slice =  cost_vol[disp,:,:]
        cur_slice = cur_slice.unsqueeze(0).unsqueeze(0)

        inputs = [leftT, cur_slice]

        test = f(*inputs)
        cost_vol[disp,:,:] = np.squeeze(test)
    return cost_vol

def createCostVol(left_im,right_im,max_disp):

    print('Creating cost-volume....')
    a_h, a_w,c = left_im.shape

    left_im = np.transpose(left_im, (2,0,1)).astype(np.uint8)
    right_im = np.transpose(right_im, (2,0,1)).astype(np.uint8)

    left_im = np.reshape(left_im, [1,c,a_h,a_w])
    right_im = np.reshape(right_im, [1,c,a_h,a_w])

    with torch.no_grad():

        left_imT = Variable(Tensor(left_im.astype(np.uint8)))
        right_imT = Variable(Tensor(right_im.astype(np.uint8)))

        left_feat = branch(left_imT)
        right_feat = branch(right_imT)

        _,f,h,w = left_feat.shape

        cost_vol = np.zeros((max_disp+1,a_h,a_w))
        cost_volT = Variable(Tensor(cost_vol))

        #0 => max_disp => one less disp!
        for disp in range(0,max_disp+1):
            if(disp == 0):
                sim_score = cos(left_feat, right_feat)
                cost_volT[disp,:,:] = torch.squeeze(sim_score)
            else:
                right_shifted = torch.cuda.FloatTensor(1,f,h,w).fill_(0)
                right_shift = torch.cuda.FloatTensor(1,f,h,disp).fill_(0)
                right_appended = torch.cat([right_shift,right_feat],3)

                _,f,h_ap,w_ap = right_appended.shape
                right_shifted[:,:,:,:] = right_appended[:,:,:,:(w_ap-disp)]
                sim_score = cos(left_feat, right_shifted)
                cost_volT[disp,:,:] = torch.squeeze(sim_score)

    print('Done')
    return cost_volT

def createCostVolRL(left_im,right_im,max_disp):

    print('Create cost-volume right-to-left...')
    a_h, a_w,c = left_im.shape

    left_im = np.transpose(left_im, (2,0,1)).astype(np.uint8)
    right_im = np.transpose(right_im, (2,0,1)).astype(np.uint8)

    left_im = np.reshape(left_im, [1,c,a_h,a_w])
    right_im = np.reshape(right_im, [1,c,a_h,a_w])

    with torch.no_grad():
        left_imT = Variable(Tensor(left_im))
        right_imT = Variable(Tensor(right_im))

        left_feat = branch(left_imT)
        right_feat = branch(right_imT)

        _,f,h,w = left_feat.shape
        cost_vol = np.zeros((max_disp+1,a_h,a_w))

        cost_volT = Variable(Tensor(cost_vol))

        for disp in range(0,max_disp+1):
            if(disp == 0):
                sim_score = cos(right_feat, left_feat)
                cost_volT[disp,:,:] = torch.squeeze(sim_score)
            else:
                left_shifted = torch.cuda.FloatTensor(1,f,h,w).fill_(0)
                left_shift = torch.cuda.FloatTensor(1,f,h,disp).fill_(0)
                left_appended = torch.cat([left_feat,left_shift],3)

                _,f,h_ap,w_ap = left_appended.shape
                left_shifted[:,:,:,:] = left_appended[:,:,:,disp:w_ap]
                sim_score = cos(right_feat, left_shifted)
                cost_volT[disp,:,:] = torch.squeeze(sim_score)

    print('Done')
    return cost_volT

def LR_Check(first_output, second_output):

    print('Check for inconsistencies...')
    h,w = first_output.shape

    line = np.array(range(0, w))
    idx_arr = np.matlib.repmat(line,h,1)

    dif = idx_arr - first_output

    first_output[np.where(dif <= 0)] = 0

    first_output = first_output.astype(np.int)
    second_output = second_output.astype(np.int)
    dif = dif.astype(np.int)

    second_arr_reordered = np.array(list(map(lambda x, y: y[x], dif, second_output)))

    dif_LR = np.abs(second_arr_reordered - first_output)
    first_output[np.where(dif_LR >= 1.1)] = 0

    first_output = first_output.astype(np.float32)
    first_output[np.where(first_output == 0.0)] = np.nan
    print('Done')
    return first_output

def FillIncons(mask, disp):

    #limit for consistent point search
    max_search = 300
    w = mask.shape[1]
    h = mask.shape[0]

    #BG
    idc = np.argwhere(np.isnan(disp))
    for curnan in range(len(idc)):
        curnanh = idc[curnan][0]
        curnanw = idc[curnan][1]
        if(mask[curnanh,curnanw] == 0):

            #whole scanline is nan => disp is 0
            if(all(np.isnan(disp[curnanh,:]))):
                #hole line set to 0!
                disp[curnanh,:] = 0.0
            #all px to the left are NaN
            if(all(np.isnan(disp[curnanh,0:curnanw]))):
                #go to the right
                curw = curnanw
                fill = 0
                while(np.isnan(disp[curnanh,curw]) and mask[curnanh,curnanw] == 0):
                    curw = curw +1
                    fill = disp[curnanh,curw]
                disp[curnanh,curnanw] = fill
            #else go left
            else:
                curw = curnanw
                fill = 0
                while(np.isnan(disp[curnanh,curw]) and mask[curnanh,curnanw] == 0):
                    curw = curw -1
                    fill = disp[curnanh,curw]
                disp[curnanh,curnanw] = fill
    #FG
    idcFG = np.argwhere(np.isnan(disp))
    for curnan in range(len(idcFG)):

        curnanh = idcFG[curnan][0]
        curnanw = idcFG[curnan][1]

        left = 0
        right = 0
        above = 0
        under = 0

        r_above = 0
        l_above = 0
        r_under = 0
        l_under = 0

        if(curnanw == 0):
            left = 0
        else:
            left = int(disp[curnanh,curnanw-1])
        counter = 0
        while(np.isnan(disp[curnanh,curnanw+counter])):
            counter = counter +1
            if((curnanw+counter) >= w or counter >= max_search):
                right = 0
                break
            right = disp[curnanh,curnanw+counter]
        counter = 0
        while(np.isnan(disp[curnanh+counter,curnanw])):
            counter = counter +1
            if((curnanh+counter) >= h or counter >= max_search):
                above = 0
                break
            above = disp[curnanh+counter,curnanw]
        if(curnanh == 0):
            under = 0
        else:
            under = disp[curnanh-1,curnanw]

        counter = 0
        while(np.isnan(disp[curnanh+counter,curnanw+counter])):
            counter = counter +1
            if((curnanh+counter) >= h or counter >= max_search):
                r_above = 0
                break
            if((curnanw+counter) >= w):
                r_above = 0
                break
            r_above = disp[curnanh+counter,curnanw+counter]

        if(curnanh == 0 or curnanw == 0):
            l_under = 0
        else:
            l_under = disp[curnanh-1,curnanw-1]

        counter = 0
        while(np.isnan(disp[curnanh+counter,curnanw-counter])):
            counter = counter +1
            if((curnanh+counter) >= h):
                l_above = 0
                break
            if((curnanw-counter) <= 0 or counter >= max_search):
                l_above = 0
                break
            l_above = disp[curnanh+counter,curnanw-counter]

        if(curnanh == 0 or curnanw >= w-1):
            r_under = 0
        else:
            r_under = disp[curnanh-1,curnanw+1]

        fill = np.median([left,right,above,under,r_above,l_above,r_under,l_under])
        disp[curnanh,curnanw] = fill
    return disp

def TestImage(fn_left, fn_right, max_disp, im_to_save, filtered = True, lr_check = True, fill = True):

    left = cv2.imread(fn_left)
    right = cv2.imread(fn_right)
    disp_map = []

    if(filtered):

        cost_vol = createCostVol(left,right,max_disp)
        cost_vol = filterCostVolMedian(cost_vol)

        cost_vol_filteredn = filterCostVolBilat(cost_vol,left)
        cost_vol_filteredn = np.squeeze(cost_vol_filteredn.cpu().data.numpy())

        disp = np.argmax(cost_vol_filteredn, axis=0)
        writePFM(im_to_save + '.pfm', disp.astype(np.float32), scale=1)
        if(lr_check):
            cost_vol_RL = createCostVolRL(left,right,max_disp)
            cost_vol_RL = filterCostVolMedian(cost_vol_RL)

            cost_vol_RL_fn = filterCostVolBilat(cost_vol_RL,right)
            cost_vol_RL_fn = np.squeeze(cost_vol_RL_fn.cpu().data.numpy())

            disp_map_RL = np.argmax(cost_vol_RL_fn, axis=0)
            disp_map = LR_Check(disp.astype(np.float32), disp_map_RL.astype(np.float32))
            writePFM(im_to_save + '_s.pfm', disp_map.astype(np.float32), scale=1)

    else:

        cost_vol = createCostVol(left,right,max_disp)
        cost_vol = np.squeeze(cost_vol.cpu().data.numpy())

        disp = np.argmax(cost_vol, axis=0)
        writePFM(im_to_save + '.pfm', disp.astype(np.float32), scale=1)

        if(lr_check):

            cost_vol_RL = createCostVolRL(left,right,max_disp)
            cost_vol_RL = np.squeeze(cost_vol_RL.cpu().data.numpy())
            disp_map_RL = np.argmax(cost_vol_RL, axis=0)
            disp_map = LR_Check(disp.astype(np.float32), disp_map_RL.astype(np.float32))
            writePFM(im_to_save + '_s.pfm', disp_map.astype(np.float32), scale=1)

    if(fill):

        disp = np.array(disp_map)
        im_disp = Image.fromarray(disp)
        im_disp = np.dstack((im_disp, im_disp, im_disp)).astype(np.uint8)

        h,w = disp.shape

        shifted = cv2.pyrMeanShiftFiltering(im_disp, 7, 7)

        gray = cv2.cvtColor(shifted, cv2.COLOR_BGR2GRAY)
        thresh = cv2.threshold(gray, 0, 1,
        cv2.THRESH_BINARY | cv2.THRESH_OTSU)[1]

        cv2.imwrite(im_to_save + 'bilat_and_med_mask.png',thresh * 255)

        disp_filled = FillIncons(thresh, disp)
        writePFM(im_to_save + '_filled.pfm',disp_filled)
    return disp_map, disp

if __name__ == "__main__":
    main()

ext: py

reason: import numpy


text:
import argparse
import numpy as np
from sklearn.metrics.pairwise import rbf_kernel, linear_kernel
from util import *

parser = argparse.ArgumentParser()
parser.add_argument("-i", "--input-dir", help="Input directory", required=True)
parser.add_argument("-b", "--batch-col", help="Batch column", required=True)
parser.add_argument("-o", "--output-path", help="Path to output file", required=True)
args = parser.parse_args()

def mmd(sample1, sample2, kernel=rbf_kernel):
    """ "the kernel we used is a sum of three Gaussian kernels with
    different scales...We chose the [sigma_i]s to be m/2,m,2m, where m
    is the median of the average distance between a point in the
    target sample to its nearest 25 neighbors, a popular practice in
    kernel-based methods."
    - https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5870543/
    """
    k_x_x = kernel(sample1, sample1)
    k_x_y = kernel(sample1, sample2)
    k_y_y = kernel(sample2, sample2)

    n = len(sample1)
    m = len(sample2)

    mean_xx = (1 / n**2) * k_x_x.sum()
    mean_xy = (1 / (m * n)) * k_x_y.sum()
    mean_yy = (1 / m**2) * k_y_y.sum()
    mmd_squared = mean_xx - 2 * mean_xy + mean_yy

    return mmd_squared ** 0.5

def mmd_multi_batch(batches):
    mmds = []
    for i in range(len(batches)):
        for j in range(i+1, len(batches)):
            mmds.append(mmd(batches[i], batches[j]))
    return np.array(mmds).mean()

def calculate_mmd(df, batch_col):
    return mmd_multi_batch(split_into_batches(df, batch_col))

# Open files
cache = DataFrameCache()

if not os.path.exists(args.output_path):
    with open(args.output_path, "w") as output_file:
        output_file.write("metric,adjuster,dataset,value\n")

# Calculate metrics
unadjusted_path = args.input_dir + "/unadjusted.csv"
unadj = cache.get_dataframe(unadjusted_path)
dataset = os.path.basename(args.input_dir)

for method in ["scaled", "combat", "confounded"]:
    adjusted_path = args.input_dir + "/" + method + ".csv"
    df = cache.get_dataframe(adjusted_path)

    print("Calculating MMD for the {} method on the {} dataset".format(method, dataset), flush=True)
    value = calculate_mmd(df, args.batch_col)

    with open(args.output_path, "a") as output_file:
        output_file.write("{},{},{},{}\n".format("MMD", method, dataset, value))


text:
from __future__ import absolute_import
from __future__ import print_function
from __future__ import division

import sys
import os
import os.path as osp

sys.path.append(osp.join(osp.dirname(osp.dirname(osp.dirname(__file__))), 'configs'))
from path_cfg import MOTSYNTH_ROOT, MOTCHA_ROOT

from torchreid.data import ImageDataset
import tqdm

import pandas as pd
import numpy as np
import json

def read_json(path):
    with open(path) as json_file:
        data = json.load(json_file)
    return data

def anns2df(anns, img_dir):
    if isinstance(anns['info'], list):
        assert len(anns['info']) == 1
        isnight = anns['info']['is_night']

    else: # Messy way to recover 'is_night' from combined annotation files
        imid2isnight = {im['id']: im['is_night'] for im in anns['images']}
        imid2frame = {im['id']: im['frame_n'] for im in anns['images']}
        isnight = {ann['id']: imid2isnight[ann['image_id']] for ann in anns['annotations']}

    rows = []
    for ann in tqdm.tqdm(anns['annotations']):
        # TODO: Get rid of this
        if ann['bbox'][2] > 1920 or ann['bbox'][3] > 1080:
            continue

        assert osp.exists(f"{img_dir}/{ann['id']}.png"), f"{img_dir}/{ann['id']}.png does not exist!!!"
        ##if not osp.exists(f"{img_dir}/{ann['id']}.png"):
        #    continue

        row={'path': f"{img_dir}/{ann['id']}.png",
            'model_id': int(ann['model_id']),
            'height': int(ann['bbox'][-1]),
            'width': int(ann['bbox'][-2]),
            'iscrowd': int(ann['iscrowd']),
            'isnight': int(isnight if not isinstance(isnight, dict) else isnight[ann['id']]),
            'vis' : (np.array(ann['keypoints'])[2::3]).mean(),
            'frame_n': int(imid2frame[ann['image_id']]),
            **{f'attr_{i}': int(attr_val) for i, attr_val in enumerate(ann['attributes'])}}
        rows.append(row)

    return  pd.DataFrame(rows)

def anns2df_motcha(anns, img_dir):
    # Build DF from anns
    rows = []
    for ann in tqdm.tqdm(anns['annotations']):
        row={'path': f"{osp.join(img_dir)}/{ann['id']}.png",
            'ped_id': int(ann['ped_id']),
            'height': int(ann['bbox'][-1]),
            'width': int(ann['bbox'][-2]),
            'iscrowd': int(ann['iscrowd']),
            'vis' : float(ann['vis']),
            'frame_n': int(ann['frame_n'])}
        rows.append(row)

    return  pd.DataFrame(rows)

def assign_ids(df, night_id=True, attr_indices = [0, 2, 3, 4, 7, 9, 10]): #attr_indices = [0, 2, 3, 4, 7, 8, 9, 10]):
    id_cols = ['model_id'] + [f'attr_{i}' for i in attr_indices if f'attr{i}' in df.columns]
    if night_id and 'isnight' in df.columns:
        id_cols += ['isnight']

    unique_ids_df = df[id_cols].drop_duplicates()
    unique_ids_df['reid_id'] = np.arange(unique_ids_df.shape[0])

    return  df.merge(unique_ids_df)

def clean_rows(df, min_vis, min_h, min_w, min_samples):
    # Filter by size and occlusion
    keep = (df['vis'] >= min_vis) & (df['height']>=min_h) & (df['width'] >= min_w) & (df['iscrowd']==0)
    clean_df = df[keep]

    # Keep only ids with at least MIN_SAMPLES appearances
    clean_df['samples_per_id'] = clean_df.groupby('reid_id')['height'].transform('count').values
    clean_df = clean_df[clean_df['samples_per_id']>=min_samples]

    return clean_df

def relabel_ids(df):
    df.rename(columns = {'reid_id': 'reid_id_old'}, inplace=True)

    # Relabel Ids from 0 to N-1
    ids_df = df[['reid_id_old']].drop_duplicates()
    ids_df['reid_id'] = np.arange(ids_df.shape[0])
    df = df.merge(ids_df)
    return df

class MOTSeqDataset(ImageDataset):
    def __init__(self,ann_file, img_dir, min_vis=0.25, min_h=50, min_w=25, min_samples=15, night_id=True, motcha=False, **kwargs):

        print("Reading json...")
        anns = read_json(ann_file)
        print("Done!")

        print("Preparing dataset...")
        if motcha:
            df = anns2df_motcha(anns, img_dir)
            df['reid_id'] = df['ped_id']

        else:
            df = anns2df(anns, img_dir)
            df = assign_ids(df, night_id=night_id, attr_indices = [0, 2, 3, 4, 7, 8, 9, 10])

        df= clean_rows(df, min_vis, min_h=min_h, min_w=min_w, min_samples=min_samples)
        df = relabel_ids(df)

        # For testing, choose one apperance randomly for every track and put in the gallery
        to_tuple_list = lambda df: list(df[['path', 'reid_id', 'cam_id']].itertuples(index=False, name=None))
        df['cam_id'] = 0
        train = to_tuple_list(df)

        df['index'] = df.index.values
        np.random.seed(0)
        query_per_id = df.groupby('reid_id')['index'].agg(lambda x: np.random.choice(list(x.unique())))
        query_df = df.loc[query_per_id.values].copy()
        gallery_df = df.drop(query_per_id).copy()
        gallery_df['cam_id'] = 1

        query=to_tuple_list(query_df)
        gallery=to_tuple_list(gallery_df)

        print("Done!")
        super(MOTSeqDataset, self).__init__(train, query, gallery, **kwargs)

def get_sequence_class(seq_name):
    """Hacky way to reuse the same class code for different sequences. We wrap the same class and
    change the __name__ attribute.
    """

    if 'MOT17' in seq_name:
        ann_file = osp.join(MOTCHA_ROOT, f'motcha_coco_annotations/{seq_name}.json')
        img_dir = osp.join(MOTCHA_ROOT, 'reid_images')
        min_samples = 5
        motcha=True

    else:
        #assert seq_name in [f'motsynth_{split}{maybe_mini}' for split in [1, 2, 3, 4, 'train', 'val'] for maybe_mini in ('_mini', '')]
        split_name = seq_name.split('motsynth_')[-1]
        ann_file = osp.join(MOTSYNTH_ROOT, 'comb_annotations', f'{split_name}.json')
        img_dir = osp.join(MOTSYNTH_ROOT, 'reid')
        min_samples = 15
        motcha=False

    class MOTSeq(MOTSeqDataset):
        def __init__(self, **kwargs):
            super(MOTSeq, self).__init__(ann_file = ann_file, img_dir=img_dir, min_samples=min_samples, motcha=motcha,**kwargs)

    MOTSeq.__name__=seq_name

    return MOTSeq

max_stars_repo_name: gritYCDA/motsynth-baselines

max_issues_repo_name: gritYCDA/motsynth-baselines

max_forks_repo_name: gritYCDA/motsynth-baselines


text:
import matplotlib.pyplot as plt
import numpy as np
import random

#加载数据
def loadDataSet():
	dataMat = [] #创建数据列表
	labelMat = [] #创建标签列表
	fr = open('testSet.txt') #打开文件
	for line in fr.readlines(): #逐行读取
		lineArr = line.strip().split() #去回车，放入列表
		dataMat.append([1.0,float(lineArr[0]),float(lineArr[1])]) #添加数据
		labelMat.append(int(lineArr[2])) #添加标签
	fr.close() #关闭文件
	return dataMat,labelMat

#sigmoid函数
def sigmoid(inX):
	return 1.0/(1+np.exp(-inX))

#梯度上升法
def gradAscent(dataMatIn,classLabels):
	dataMatrix = np.mat(dataMatIn) #转换成numpy的mat
	labelMat = np.mat(classLabels).transpose() #转换成numpy的mat,并进行转置
	m,n = np.shape(dataMatrix) #返回dataMatrix的大小，m为行数，n为列数
	alpha = 0.01 #移动步长，用来控制更新的幅度
	maxCycles = 500 #最大迭代次数
	weights = np.ones((n,1)) #每个回归系数初始化为1
	weights_array = np.array([])
	for k in range(maxCycles): #梯度上升矢量化公式
		h = sigmoid(dataMatrix*weights)
		error = labelMat - h
		weights = weights + alpha*dataMatrix.transpose()*error
		weights_array = np.append(weights_array,weights)
	weights_array = weights_array.reshape(maxCycles,n)
	return weights.getA(),weights_array #将矩阵转换为数组,返回权重数组 weights.getA()求得的权重数组(最优参数), weights_array 每次更新的回归系数

#随机梯度上升法
def stocGradAscent(dataMatrix,classLabels,numIter=150):
	m,n = np.shape(dataMatrix) #返回dataMatrix的大小，m为行数，n为列数
	weights = np.ones(n)  #参数初始化
	weights_array = np.array([])  #存储每次更新的回归系数
	for j in range(numIter):
		dataIndex = list(range(m))
		for i in range(m):
			alpha = 4/(1.0+j+i)+0.01   #降低alpha的大小，每次减小1/(j+i)
			randIndex = int(random.uniform(0,len(dataIndex))) #随机选取样本
			h = sigmoid(sum(dataMatrix[randIndex]*weights)) #选择随机选取的一个样本，计算h
			error = classLabels[randIndex]-h #计算误差
			weights = weights+alpha*error*dataMatrix[randIndex]  #更新回归系数
			weights_array = np.append(weights_array,weights,axis=0)  #添加回归系数到数组中
			del(dataIndex[randIndex]) #删除已经使用的样本
	weights_array = weights_array.reshape(numIter*m,n) #改变纬度
	return weights,weights_array   #weights - 求得的回归系数数组(最优参数) weights_array - 每次更新的回归系数

#绘制回归系数与迭代次数的关系
def plotWeights(weights_array1,weights_array2):
	fig,axs = plt.subplots(nrows=3,ncols=2,sharex=False,sharey=False,figsize=(13,8))   #当nrow=3,nclos=2时,代表fig画布被分为六个区域,axs[0][0]表示第一行第一列
	x1 = np.arange(0,len(weights_array1),1)
	#绘制w0与迭代次数的关系
	axs[0][0].plot(x1,weights_array1[:,0])
	axs0_title_text = axs[0][0].set_title(u'Stochastic Gradient Ascent Algorithm')
	axs0_ylabel_text = axs[0][0].set_ylabel(u'W0')
	plt.setp(axs0_title_text,size=18,color='black')
	plt.setp(axs0_ylabel_text,size=18,color='black')
	#绘制w1与迭代次数的关系
	axs[1][0].plot(x1,weights_array1[:,1])
	axs1_ylabel_text = axs[1][0].set_ylabel(u'W1')
	plt.setp(axs1_ylabel_text,size=18,color='black')
	#绘制w2与迭代次数的关系
	axs[2][0].plot(x1,weights_array1[:,2])
	axs2_xlabel_text = axs[2][0].set_xlabel(u'Iteration Number')
	axs2_ylabel_text = axs[2][0].set_ylabel(u'W2')
	plt.setp(axs2_xlabel_text,size=18,color='black')
	plt.setp(axs2_ylabel_text,size=18,color='black')

	x2 = np.arange(0,len(weights_array2),1)
	#绘制w0与迭代次数的关系
	axs[0][1].plot(x2,weights_array2[:,0])
	axs0_title_text = axs[0][1].set_title(u'Gradient Ascent Algorithm')
	axs0_ylabel_text = axs[0][1].set_ylabel(u'W0')
	plt.setp(axs0_title_text,size=18,color='black')
	plt.setp(axs0_ylabel_text,size=18,color='black')
	#绘制w1与迭代次数的关系
	axs[1][1].plot(x2,weights_array2[:,1])
	axs1_ylabel_text = axs[1][1].set_ylabel(u'W1')
	plt.setp(axs1_ylabel_text,size=18,color='black')
	#绘制w2与迭代次数的关系
	axs[2][1].plot(x2,weights_array2[:,2])
	axs2_xlabel_text = axs[2][1].set_xlabel(u'Iteration Number')
	axs2_ylabel_text = axs[2][1].set_ylabel(u'W2')
	plt.setp(axs2_xlabel_text,size=18,color='black')
	plt.setp(axs2_ylabel_text,size=18,color='black')

	plt.show()

if __name__=='__main__':
	dataMat,labelMat = loadDataSet()
	weights1,weights_array1 = stocGradAscent(np.array(dataMat),labelMat)
	weights2,weights_array2 = gradAscent(dataMat,labelMat)
	plotWeights(weights_array1,weights_array2)

max_stars_repo_name: BettyMei/ML-python

max_issues_repo_name: BettyMei/ML-python

max_forks_repo_name: BettyMei/ML-python


text:
import os
import numpy as np
import matplotlib.pyplot as plt
import tensorflow as tf
import seq2seq

def plot(data):
    framerate = 60
    points = 2000
    skip = int(len(data) / points)
    duration = len(data) / float(framerate)
    data_skipped = data.transpose(1, 0, 2)[:, ::skip, :]
    plt.figure(figsize=(16, 10), facecolor='white')
    for i, joint in enumerate(data_skipped):
        plt.gca().set_prop_cycle(None)
        plt.plot(np.linspace(0, duration, len(joint)), i + joint, lw=1)
    plt.ylim([0, data_skipped.shape[0]])
    plt.xlim([0, duration])
    plt.ylabel('Joint')
    plt.xlabel('Seconds')
    plt.show()

def fixed_generator(data, hop_length, sequence_length, batch_size=50):
    idxs = np.arange(0, len(data) - sequence_length, hop_length)
    n_batches = len(idxs)
    for batch_i in range(n_batches):
        this_idxs = idxs[batch_i * batch_size:(batch_i + 1) * batch_size]
        source = [data[i:i + sequence_length, :] for i in this_idxs]
        target = [
            data[i + sequence_length:i + sequence_length * 2, :]
            for i in this_idxs
        ]
        yield np.array(
            source, dtype=np.float32), np.array(
                target, dtype=np.float32)

def batch_generator(data, sequence_length, batch_size=50):
    idxs = np.random.permutation(np.arange(len(data) - sequence_length * 2))
    n_batches = len(idxs) // (batch_size * sequence_length)
    for batch_i in range(n_batches):
        this_idxs = idxs[batch_i * batch_size:(batch_i + 1) * batch_size]
        source = [data[i:i + sequence_length, :] for i in this_idxs]
        target = [
            data[i + sequence_length:i + sequence_length * 2, :]
            for i in this_idxs
        ]
        yield np.array(
            source, dtype=np.float32), np.array(
                target, dtype=np.float32)

def train(data,
          data_mean,
          data_std,
          n_epochs=1000,
          batch_size=100,
          sequence_length=240,
          ckpt_path='./',
          model_name='seq2seq.ckpt',
          restore_name=None,
          **kwargs):
    config = tf.ConfigProto()
    config.gpu_options.allow_growth = True
    sess = tf.Session(config=config)

    net = seq2seq.create_model(
        batch_size=batch_size, sequence_length=sequence_length, **kwargs)

    learning_rate = tf.placeholder(tf.float32, name='learning_rate')
    opt = tf.train.AdamOptimizer(learning_rate=learning_rate).minimize(
        net['loss'])
    init_op = tf.group(tf.global_variables_initializer(),
                       tf.local_variables_initializer())
    sess.run(init_op)
    saver = tf.train.Saver()
    if restore_name:
        saver.restore(sess, restore_name)

    current_learning_rate = 0.0001
    for epoch_i in range(n_epochs):
        total, total_mse, total_mdn = 0.0, 0.0, 0.0
        for it_i, (source, target) in enumerate(
                batch_generator(
                    data,
                    sequence_length=sequence_length,
                    batch_size=batch_size)):
            mse_loss, mdn_loss, _ = sess.run(
                [net['mse_loss'], net['mdn_loss'], opt],
                feed_dict={
                    learning_rate: current_learning_rate,
                    net['keep_prob']: 0.8,
                    net['source']: source,
                    net['target']: target
                })
            total += mse_loss + mdn_loss
            total_mdn += mdn_loss
            total_mse += mse_loss
            print(
                '{}: mdn: {}, mse: {}, total: {}'.format(
                    it_i, mdn_loss, mse_loss, mdn_loss + mse_loss),
                end='\r')
        current_learning_rate = max(0.0001, current_learning_rate * 0.99)
        print('iteration: {}, learning rate: {}'.format(it_i,
                                                        current_learning_rate))
        print('\n-- epoch {}: mdn: {}, mse: {}, total: {} --\n'.format(
            epoch_i, total_mdn / (it_i + 1), total_mse / (it_i + 1),
            total / (it_i + 1)))
        saver.save(
            sess, os.path.join(ckpt_path, model_name), global_step=epoch_i)

    sess.close()

def infer(source,
          target,
          data_mean,
          data_std,
          batch_size,
          sequence_length,
          ckpt_path='./',
          model_name='seq2seq.ckpt',
          **kwargs):
    config = tf.ConfigProto()
    config.gpu_options.allow_growth = True
    with tf.Graph().as_default() as g, tf.Session(config=config) as sess:
        net = seq2seq.create_model(
            batch_size=batch_size, sequence_length=sequence_length, **kwargs)

        init_op = tf.group(tf.global_variables_initializer(),
                           tf.local_variables_initializer())
        sess.run(init_op)
        saver = tf.train.Saver()
        saver.restore(sess, os.path.join(ckpt_path, model_name))
        recon, enc = sess.run(
            [net['decoding'], net['encoding']],
            feed_dict={
                net['source']: source,
                net['keep_prob']: 1.0
            })
        src = (source * data_std) + data_mean
        tgt = (target * data_std) + data_mean
        res = np.minimum(1.0, np.maximum(0.0, (recon[1] * data_std) + data_mean))
        fig, axs = plt.subplots(1, 3, sharey=True)
        axs[0].plot(src.reshape(-1, src.shape[-1]))
        axs[0].set_title('Source')
        axs[1].plot(tgt.reshape(-1, tgt.shape[-1]))
        axs[1].set_title('Target (Original)')
        axs[2].plot(res.reshape(-1, res.shape[-1]))
        axs[2].set_title('Target (Synthesis Sampling)')
        return {
            'source': src,
            'target': tgt,
            'encoding': enc,
            'prediction': res
        }

# def feedback(source,
#              target,
#              data_mean,
#              data_std,
#              batch_size,
#              sequence_length,
#              ckpt_path='./',
#              model_name='seq2seq.ckpt',
#              **kwargs):
#     config = tf.ConfigProto()
#     config.gpu_options.allow_growth = True
#     with tf.Graph().as_default() as g, tf.Session(config=config) as sess:
#         net = seq2seq.create_model(
#             batch_size=batch_size, sequence_length=sequence_length, **kwargs)
#
#         init_op = tf.group(tf.global_variables_initializer(),
#                            tf.local_variables_initializer())
#         sess.run(init_op)
#         saver = tf.train.Saver()
#         saver.restore(sess, os.path.join(ckpt_path, model_name))
#         recon, enc = sess.run(
#             [net['decoding'], net['encoding']],
#             feed_dict={
#                 net['source']: source,
#                 net['keep_prob']: 1.0
#             })
#         src = (source * data_std) + data_mean
#         tgt = (target * data_std) + data_mean
#         res = np.minimum(1.0, np.maximum(0.0, (recon[1] * data_std) + data_mean))
#         fig, axs = plt.subplots(1, 3, sharey=True)
#         axs[0].plot(src.reshape(-1, src.shape[-1]))
#         axs[0].set_title('Source')
#         axs[1].plot(tgt.reshape(-1, tgt.shape[-1]))
#         axs[1].set_title('Target (Original)')
#         axs[2].plot(res.reshape(-1, res.shape[-1]))
#         axs[2].set_title('Target (Synthesis Sampling)')
#         np.save('source.npy', src)
#         np.save('target.npy', tgt)
#         np.save('encoding.npy', enc)
#         np.save('prediction.npy', res)
#         return {
#             'source': src,
#             'target': tgt,
#             'encoding': enc,
#             'prediction': res
#         }


text:
"""
Model use utilities.
"""
import os

import numpy as np
import tensorflow as tf
import tensorflow_hub as hub

MODEL_URL = "https://tfhub.dev/google/humpback_whale/1"

# Model saved locally here
LOCAL_MODEL = "google/humpback_whale/1"

# We apply the model to get a 1-sec score resolution
# (Note that the model expects input signals sampled at 10 kHz.)
CONTEXT_STEP_SAMPLES = tf.cast(10_000, tf.int64)

class ModelHelper:
    """
    Helps with loading and applying the Google model.
    """

    def __init__(self):
        self.model = None
        self.score_fn = None

    def load_model(self) -> None:
        """
        Loads the model, initially by downloading it from its original place,
        then from a local copy in subsequent calls.
        The model and score function are maintained internally.
        """

        if os.path.isdir(LOCAL_MODEL):
            print(f"\n==> Loading model from {LOCAL_MODEL}")
            model = hub.load(LOCAL_MODEL)
        else:
            print(f"\n==> Loading model from {MODEL_URL}")
            model = hub.load(MODEL_URL)
            print(f"    Saving model to {LOCAL_MODEL} ...")
            tf.saved_model.save(
                model,
                LOCAL_MODEL,
                signatures={"score": model.score, "metadata": model.metadata},
            )
            print(f"    Model saved to {LOCAL_MODEL}")

        self.model = model

        self.score_fn = model.signatures["score"]

        metadata_fn = model.signatures["metadata"]
        metadata = metadata_fn()
        print("metadata:")
        for key, val in metadata.items():
            print(f"  {key}: {val}")

    def apply_model(self, psound: np.ndarray) -> np.ndarray:
        """
        Applies the model on given audio segment.

        :param psound:  The input signal, assumed sampled at 10 kHz.
        :return:        Array of corresponding scores at 1-sec resolution.
        """
        waveform1 = np.expand_dims(psound, axis=1)
        waveform_exp = tf.expand_dims(waveform1, 0)  # makes a batch of size 1

        psound_scores = self.score_fn(
            waveform=waveform_exp, context_step_samples=CONTEXT_STEP_SAMPLES
        )
        return psound_scores["scores"].numpy()[0]

max_stars_repo_name: mbari-org/humpback-whale-song-detection

max_issues_repo_name: mbari-org/humpback-whale-song-detection

max_forks_repo_name: mbari-org/humpback-whale-song-detection


text:
#!python3
## create and write depth file for HYCOM
## include utilities function
import sys
import os
from os.path import join
### might not need next line in older or newer version of proplot
os.environ['PROJ_LIB'] = '/Users/abozec/opt/anaconda3/share/proj'
import proplot as plot
import numpy as np

iodir='/Users/abozec/Documents/GitHub/BB86_PACKAGE/PYTHON/'
sys.path.append(iodir)
from hycom.io import write_hycom_grid

# define path
io=join(iodir+'../topo/')
file_grid_new='regional.grid.BB86.python'

## size of the domain
idm = 101 ; jdm = 101

## longitude/latitude starting point + resolution  (in degrees) (NB: those
## variables are  not used in HYCOM)
ini_lon = 0. ; ini_lat = 0.
res = 0.20

## scale dx and dy (in m) (used in HYCOM)
dx = 20e3 ; dy = dx

## grid type
mapflg=0  ##  uniform or mercator

######### END of the USER inputs #############
## missing value in HYCOM
vmiss = 2.**100

## Get the p-point the grid
plon = np.zeros([jdm, idm])
plat = np.zeros([jdm, idm])

## Longitude
plon[:, 0] = ini_lon
for i in np.arange(idm-1)+1:
   plon[:,i] = plon[:, i-1] + res

 ## Latitude
plat[0, :] = ini_lat
for j in np.arange(jdm-1)+1:
   plat[j, :] = plat[j-1, :] + res

print('p-points grid OK')

## Declaration of the grid tabs
qlon = np.zeros([jdm,idm])
qlat = np.zeros([jdm,idm])
ulon = np.zeros([jdm,idm])
ulat = np.zeros([jdm,idm])
vlon = np.zeros([jdm,idm])
vlat = np.zeros([jdm,idm])

pang = np.zeros([jdm,idm])

pscx = np.zeros([jdm,idm])
pscy = np.zeros([jdm,idm])

qscx = np.zeros([jdm,idm])
qscy = np.zeros([jdm,idm])

uscx = np.zeros([jdm,idm])
uscy = np.zeros([jdm,idm])

vscx = np.zeros([jdm,idm])
vscy = np.zeros([jdm,idm])

cori = np.zeros([jdm,idm])
pasp = np.zeros([jdm,idm])

## Longitude/Latitude for each point (Q-points stay at the same location)
## longitude
vlon = plon
for i in np.arange(idm-1)+1:
   qlon[:, i] = 0.5*(plon[:,i]+plon[:, i-1])

diff = plon[0, 3]-plon[0, 2]
qlon[:, 0] = plon[:, 0]-0.5*diff
ulon = qlon

## latitude
ulat = plat
for j in np.arange(jdm-1)+1:
   qlat[j, :] = 0.5*(plat[j, :]+plat[j-1, :])

diff = plat[3, 0]-plat[2, 0]
qlat[0, :] = plat[0, :]-0.5*diff
vlat = qlat

## simplified grid with prescribed dx and dy
pscx[:,:] = dx ; pscy[:,:] = dy
uscx[:,:] = dx ; uscy[:,:] = dy
vscx[:,:] = dx ; vscy[:,:] = dy
qscx[:,:] = dx ; qscy[:,:] = dy

## Coriolis
beta = 2.e-11
for  j in np.arange(jdm):
   for  i in np.arange(idm):
      cori[j, i]=.93e-4+(j-int((jdm-2)/2))*dx*beta

## pang ( p-angle  for rotated grid)
## here uniform or mercator so pang = 0.

## pasp= paspect: pscx/pscy
pasp=pscx/pscy

##  Writing new grid file
write_hycom_grid(join(io+file_grid_new+'.a'), plon, plat, qlon,  qlat,  ulon, ulat, \
 vlon, vlat, pang,pscx, pscy, qscx,  qscy, uscx, uscy,  vscx, vscy, cori, pasp,mapflg)

print('Writing grid file done ')


text:
import torch
import torch.nn as nn
from torch.nn import functional as F
import numpy as np
import sys

class BasicBlock(nn.Module):
	"""
	Basic block is composed of 2 CNN layers with residual connection.
	Each CNN layer is followed by batchnorm layer and swish activation
	function.
	Args:
		in_channel: number of input channels
		out_channel: number of output channels
		k: (default = 1) kernel size
	"""
	def __init__(self, in_channel, out_channel, k=1):
		super(BasicBlock, self).__init__()
		self.conv1 = nn.Conv2d(
			in_channel,
			out_channel,
			kernel_size=k,
			padding=(0, 0),
			stride=(1, 1))
		self.bn1 = nn.BatchNorm2d(out_channel)

		self.conv2 = nn.Conv2d(
			out_channel,
			out_channel,
			kernel_size=1,
			padding=(0, 0),
			stride=(1, 1))
		self.bn2 = nn.BatchNorm2d(out_channel)

		self.shortcut = nn.Sequential()
		# if in_channel != out_channel:
		self.shortcut.add_module(
			'conv',
			nn.Conv2d(
				in_channel,
				out_channel,
				kernel_size=k,
				padding=(0,0),
				stride=(1,1)))
		self.shortcut.add_module('bn', nn.BatchNorm2d(out_channel))

	def swish(self,x):
		"""
		We use swish in spatio-temporal encoding/decoding. We tried with
		other activation functions such as ReLU and LeakyReLU. But we
		achieved the best performance with swish activation function.
		Args:
			X: tensor: (batch_size, ...)
		Return:
			_: tensor: (batch, ...): applies swish
			activation to input tensor and returns
		"""
		return x*torch.sigmoid(x)

	def forward(self, x):
		y = self.swish(self.bn1(self.conv1(x)))
		y = self.swish(self.bn2(self.conv2(y)))
		y = y + self.shortcut(x)
		y = self.swish(y)
		return y

class Model(nn.Module):
	"""
	Model for our proposed GMVAE model for human motion generation. It is composed
	of 3 encoder modules and 3 decoder modules.
	Args:
		num_class: number of action categories
		latent_dim: bottleneck dimension of the autoencoder
		components: number of Gaussian components used for sampling. We keep the
		number of components same as the number of classes
	"""
	def __init__(self, num_class, latent_dim, components=120):
		super(Model, self).__init__()
		self.latent_dim = latent_dim

		# encoder
		self.encoder1 = BasicBlock(1,1,k=5)
		self.encode_t = BasicBlock(64, 32)
		self.encode_t1 = BasicBlock(32, 8)
		self.encode_t2 = BasicBlock(8, 4)

		# decoder
		self.conv1 = BasicBlock(1,1)
		self.decode_t = BasicBlock(4,8)
		self.decode_t1 = BasicBlock(8,32)
		self.decode_t2 = BasicBlock(32,64)
		self.decoder = nn.Linear(88, 144*2)

		# Reparameterization
		self.mean = nn.Linear(352, self.latent_dim)
		self.logvar = nn.Linear(352, self.latent_dim)
		self.hidden2latent = nn.Linear(352+20+num_class+6+4, 352*2)
		self.latent2hidden = nn.Linear(self.latent_dim+num_class+6, 352)

		# root trajectory
		self.root1 = nn.Conv1d(4,8,5)
		self.root2 = nn.Conv1d(8,16,5)
		self.root3 = nn.Conv1d(16,32,5)
		self.root4 = nn.Conv1d(32,64,5)
		self.root5 = nn.Linear(72,3*2)

		# root trajectory encoder
		self.r_encoder0 = nn.Conv1d(64, 32, 1)
		self.r_encoder1 = nn.Conv1d(32, 16, 1)
		self.r_encoder2 = nn.Conv1d(16, 8, 1)
		self.r_encoder3 = nn.Conv1d(8, 4, 1)
		self.r_encoder4 = nn.Linear(24, 20)
		# self.r_encoder4 = nn.Linear(32, 16, 5)

		# sequence length encoder
		self.seq_encoder0 = nn.Conv1d(64,32,1)
		self.seq_encoder1 = nn.Conv1d(32,16,1)
		self.seq_encoder2 = nn.Conv1d(16,8,1)
		self.seq_encoder3 = nn.Conv1d(8,4,1)
		self.seq_encoder4 = nn.Linear(4, 4)

		# sequence length decoder
		# decoder
		self.seq_decoder1 = nn.Linear(latent_dim, 4)
		self.seq_decoder2 = nn.Conv1d(4,8,1)
		self.seq_decoder3 = nn.Conv1d(8,16,1)
		self.seq_decoder4 = nn.Conv1d(16,32,1)
		self.seq_decoder5 = nn.Conv1d(32,64,1)
		self.relu = nn.ReLU()
		self.sigmoid = nn.Sigmoid()

		# Gausiam mixture parameters
		self.components = components

		# mixture of Gaussian parameters
		self.z_pre = torch.nn.Parameter(torch.randn(1, 2 * self.components, self.latent_dim)
                                        / np.sqrt(self.components * self.latent_dim))

		# Uniform weighting
		self.pi = torch.nn.Parameter(torch.ones(components) / components, requires_grad=False)

	def swish(self,x):
		"""
		We use swish in spatio-temporal encoding/decoding. We tried with
		other activation functions such as ReLU and LeakyReLU. But we
		achieved the best performance with swish activation function.
		Args:
			X: tensor: (batch_size, ...)
		Return:
			_: tensor: (batch, ...): applies swish
			activation to input tensor and returns
		"""
		return x*torch.sigmoid(x)

	def encoder_net(self, X):
		"""
		Encoder first downsamples the input motion in the spatial dimension
		and then downsamples in the temporal dimension and returns spatio-
		temporal feature.
		Args:
			X: tensor: (batch_size, 64, 48, 6): input motion of 2 persons. 24
			joints for each persons so total 48 joints.
		Return:
			x: tensor: (batch_size, 4, ...): spatio-temporal feature
		"""
		N,T,J = X.shape
		# pose encoding
		x = X.reshape((N*T,1,48,6))
		x = self.encoder1(x)

		# temporal encoding
		x = x.reshape((N,T,44,2))
		x = self.encode_t(x)
		x = self.encode_t1(x)
		x = self.encode_t2(x)
		x = x.reshape((N,4,-1))

		return x

	def decoder_net(self, X):
		"""
		The deocder is opposit of the encoder. It takes the vector sampled
		from a mixture of gaussian parameter conditioned by class label on-
		hot vector and viewpoint vector, upsamples it in the temporal dimension
		first and then upsamples it in the spatial dimension.
		Args:
			X: tensor: (batch_size, 4, ...): sampled vector conditionied on class
			label and viewpoint
		Return:
			x: tensor: (batch_size, 64, 48, 6): generated human motion
		"""
		N,T,J = X.shape
		# temporal decoding
		x = X.reshape((N,T,44,2))
		x = self.decode_t(x)
		x = self.decode_t1(x)
		x = self.decode_t2(x)

		# pose decoding
		x = x.reshape((N*64,1,44,2))
		x = self.conv1(x)
		x = x.reshape((N,64, -1))
		x = self.decoder(x)

		return x

	def root_traj(self, z):
		"""
		This function calculate the root trajectory for 2-person inteaction
		classes. generates the displacement of the second person's root from the
		first person's root.
		Args:
			z: tensor: (batch_size, 4, ...): sampled vector conditionied on class
			label and viewpoint
		Return:
			z: tensor: (batch_size, 64, 3): displacement
		"""
		z = self.swish(self.root1(z))
		z = self.swish(self.root2(z))
		z = self.swish(self.root3(z))
		z = self.swish(self.root4(z))
		z = self.root5(z)

		return z

	def root_traj_encoder(self, root):
		root = root.float()
		z = self.swish(self.r_encoder0(root))
		z = self.swish(self.r_encoder1(z))
		z = self.swish(self.r_encoder2(z))
		z = self.swish(self.r_encoder3(z))
		z = z.reshape((z.shape[0], -1))
		z = self.r_encoder4(z)

		return z

	def seq_encoder(self, x):
		N,T,_ = x.shape
		z = self.relu(self.seq_encoder0(x))
		z = self.relu(self.seq_encoder1(z))
		z = self.relu(self.seq_encoder2(z))
		z = self.relu(self.seq_encoder3(z))
		z = z.reshape((N, -1))
		z = self.relu(self.seq_encoder4(z))
		return z

	def seq_decoder(self, z):
		N,_ = z.shape
		z = self.relu(self.seq_decoder1(z))
		z = z.unsqueeze(-1)
		z = self.relu(self.seq_decoder2(z))
		z = self.relu(self.seq_decoder3(z))
		z = self.relu(self.seq_decoder4(z))
		z = self.seq_decoder5(z)
		z = self.sigmoid(torch.cumsum(z.squeeze(2), dim=1).unsqueeze(2)) # monotonically increasing sequence

		return z

	def reparameterization(self,mean,logvar):
		"""
		this function sample from a unimodal Gaussian distribution in case of
		vanilla VAE.
		Args:
			mean: tensor: (batch_size, ...): mean
			logvar: tensor: (batch_size, ...): var
		Return:
			_: tensor: (batch_size, ...): sampled latent vector
		"""
		std = torch.exp(0.5*logvar)
		eps = torch.randn_like(std).to(logvar.device)
		return eps*std + mean

	def forward(self, x, y, rot, root, seq):
		x = self.encoder_net(x.float())
		seq = self.seq_encoder(seq)
		root_encoding = self.root_traj_encoder(root)
		# reparameterization
		z = x.reshape((x.shape[0], -1))
		z = torch.cat((z,root_encoding, seq, y.float(), rot.float()), dim=1)
		z = self.hidden2latent(z)
		# mean = self.mean(z)
		# logvar = self.logvar(z)
		mean, var = self.gaussian_parameters(z, dim=1)
		# print(mean.shape, var.shape)

		# Gaussian mixture

		prior = self.gaussian_parameters(self.z_pre, dim=1)

		z = self.sample_gaussian(mean, var)

		# terms for KL divergence
		log_q_phi = self.log_normal(z, mean, var)
		log_p_theta = self.log_normal_mixture(z, prior[0], prior[1])
		kld = torch.mean(log_q_phi - log_p_theta)

		# z = self.reparameterization(mean, logvar)
		z = torch.cat((z,y.float(), rot.float()), dim=1)
		z = self.latent2hidden(z)
		seq_pred = self.seq_decoder(z)
		z = z.reshape((z.shape[0], 4, -1))

		x = self.decoder_net(z)
		root = self.root_traj(z)
		return x, kld, root, seq_pred

	def sample_gaussian(self, m, v):
		"""
		Element-wise application reparameterization trick to sample from Gaussian
		Args:
			m: tensor: (batch, ...): Mean
			v: tensor: (batch, ...): Variance
		Return:
			z: tensor: (batch, ...): Samples
		"""
		sample = torch.randn(m.shape).to(m.device)

		z = m + (v**0.5)*sample
		return z

	def gaussian_parameters(self, h, dim=-1):
		"""
		Converts generic real-valued representations into mean and variance
		parameters of a Gaussian distribution
		Args:
			h: tensor: (batch, ..., dim, ...): Arbitrary tensor
			dim: int: (): Dimension along which to split the tensor for mean and
				variance
		Returns:z
			m: tensor: (batch, ..., dim / 2, ...): Mean
			v: tensor: (batch, ..., dim / 2, ...): Variance
		"""
		m, h = torch.split(h, h.size(dim) // 2, dim=dim)
		v = F.softplus(h) + 1e-8
		return m, v

	def log_normal(self, x, m, v):
		"""
		Computes the elem-wise log probability of a Gaussian and then sum over the
		last dim. Basically we're assuming all dims are batch dims except for the
		last dim.
		Args:
			x: tensor: (batch, ..., dim): Observation
			m: tensor: (batch, ..., dim): Mean
			v: tensor: (batch, ..., dim): Variance
		Return:
			kl: tensor: (batch1, batch2, ...): log probability of each sample. Note
				that the summation dimension (dim=-1) is not kept
		"""

		const = -0.5*x.size(-1)*torch.log(2*torch.tensor(np.pi))
		log_det = -0.5*torch.sum(torch.log(v), dim = -1)
		log_exp = -0.5*torch.sum( (x - m)**2/v, dim = -1)
		log_prob = const + log_det + log_exp

		return log_prob

	def log_normal_mixture(self, z, m, v):
		"""
		Computes log probability of a uniformly-weighted Gaussian mixture.
		Args:
			z: tensor: (batch, dim): Observations
			m: tensor: (batch, mix, dim): Mixture means
			v: tensor: (batch, mix, dim): Mixture variances
		Return:
			log_prob: tensor: (batch,): log probability of each sample
		"""
		z = z.unsqueeze(1)
		log_probs = self.log_normal(z, m, v)
		log_prob = self.log_mean_exp(log_probs, 1)

		return log_prob

	def log_mean_exp(self, x, dim):
		"""
		Compute the log(mean(exp(x), dim)) in a numerically stable manner
		Args:
			x: tensor: (...): Arbitrary tensor
			dim: int: (): Dimension along which mean is computed
		Return:
			_: tensor: (...): log(mean(exp(x), dim))
		"""
		return self.log_sum_exp(x, dim) - np.log(x.size(dim))

	def log_sum_exp(self, x, dim=0):
		"""
		Compute the log(sum(exp(x), dim)) in a numerically stable manner
		Args:
			x: tensor: (...): Arbitrary tensor
			dim: int: (): Dimension along which sum is computed
		Return:
			_: tensor: (...): log(sum(exp(x), dim))
		"""
		max_x = torch.max(x, dim)[0]
		new_x = x - max_x.unsqueeze(dim).expand_as(x)
		return max_x + (new_x.exp().sum(dim)).log()

max_stars_repo_name: skelemoa/mugl

max_issues_repo_name: skelemoa/mugl

max_forks_repo_name: skelemoa/mugl


text:
""" Train
"""
from __future__ import absolute_import
from __future__ import division
from __future__ import print_function

import numpy as np
import tensorflow as tf

MAX_GEN = 40
_END_ID = 3

def reverse_decode(vector, rev_vocab):
    tokens = [rev_vocab[i] for i in vector.tolist()]
    return " ".join(tokens)

# TODO: only works for one example
# TODO: tensorflow native op
def greedy_argmax(vector, step):
    counter = 0
    while counter < MAX_GEN:
        probs = step(vector)
        last_vector = probs[-1, :]
        last_id = np.argmax(last_vector)
        vector = np.concatenate([vector, [last_id]])
        if last_id == _END_ID:
            break
        counter += 1
    return vector

def beam_search():
    pass


text:
# MIT License
#
# Copyright (c) 2019 Tuomas Halvari, Juha Harviainen, Juha Mylläri, Antti Röyskö, Juuso Silvennoinen
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

import numpy as np

from dpemu.nodes import Array, Series, TupleSeries, Tuple
from dpemu.filters.common import Missing
from dpemu.filters.time_series import SensorDrift

def test_array_works_with_regular_arrays():
    a = [0.]
    x_node = Array()
    x_node.addfilter(Missing("prob", "m_val"))
    params = {"prob": 1., "m_val": np.nan}
    out = x_node.generate_error(a, params)
    assert np.isnan(out[0])

def test_series_and_array_work_with_regular_arrays():
    a = [0.]
    x_node = Array()
    x_node.addfilter(Missing("prob", "m_val"))
    series_node = Series(x_node)
    params = {"prob": 1., "m_val": np.nan}
    out = series_node.generate_error(a, params)
    assert np.isnan(out[0])

def test_array_works_with_numpy_arrays():
    a = np.array([0.])
    x_node = Array()
    x_node.addfilter(Missing("prob", "m_val"))
    params = {"prob": 1., "m_val": np.nan}
    out = x_node.generate_error(a, params)
    assert np.isnan(out[0])

def test_tuple_series_works_with_numpy_arrays():
    a = np.array([[1, 2, 3]])
    b = np.array([[1, 1, 1]])
    data = (a, b)
    x_node = Array()
    x_node.addfilter(SensorDrift("a"))
    y_node = Array()
    root_node = TupleSeries([x_node, y_node])
    res = root_node.generate_error(data, {'a': 1})
    assert np.array_equal(res[0], np.array([[2, 4, 6]])) and np.array_equal(res[1], np.array([[1, 1, 1]]))

def test_tuple_node_works():
    data = np.array([(0., 1.), (2., 3.)])
    x_node = Tuple()
    x_node.addfilter(Missing("prob", "m_val"))
    series_node = Series(x_node)
    data = series_node.generate_error(data, {'prob': 1, 'm_val': np.nan})
    assert np.isnan(data[0][0]) and np.isnan(data[0][1]) and np.isnan(data[1][0]) and np.isnan(data[1][1])

def test_exception_raised_when_param_missing():
    data = np.random.rand(5)
    x_node = Array()
    x_node.addfilter(Missing("prob", "m_val"))
    try:
        x_node.generate_error(data, {"probb": .5, "m_val": np.nan})
    except Exception as e:
        assert "prob" in str(e)

max_forks_repo_name: thalvari/dpEmu-AutoML