text:
import os
import du

trial_name = os.path.basename(__file__)[:-3]
with du.trial.run_trial(trial_name=trial_name) as trial:

    NUM_EPOCHS = 25
    BATCH_SIZE = 500
    NUM_UNITS = 512

    import numpy as np
    import tensorflow as tf
    import tfu

    train, valid, test = du.tasks.image_tasks.mnist(x_dtype="float32",
                                                    y_dtype="int64")

    x = tf.placeholder(tf.float32, shape=[None, 28, 28, 1])
    y = tf.placeholder(tf.int64, shape=[None])

    tfu.add_hook(tfu.hooks.reuse_variables(variable_scope="valid",
                                           replace=[("valid", "train")]))
    tfu.add_hook(tfu.inits.set_weight_init(tfu.inits.msr_normal))
    tfu.add_hook(tfu.inits.scale_weight_inits(np.sqrt(2)))

    epoch_size = len(train["x"])
    batches_per_epoch = np.ceil(epoch_size / float(BATCH_SIZE))
    expected_count = NUM_EPOCHS * batches_per_epoch
    tfu.counter.make_default_counter(expected_count=expected_count)

    def model():
        h = tfu.flatten(x, 2)
        with tfu.variable_scope("fc1"):
            h = tfu.affine(h, NUM_UNITS)
            h = tf.nn.relu(h)
            h = tfu.dropout(h, keep_prob=0.5)
        with tfu.variable_scope("fc2"):
            h = tfu.affine(h, NUM_UNITS)
            h = tf.nn.relu(h)
            h = tfu.dropout(h, keep_prob=0.5)
        h = tfu.affine(h, 10, "logit")

        cross_entropy = tf.reduce_mean(
            tfu.softmax_cross_entropy_with_logits(h, y))
        accuracy = tf.reduce_mean(tfu.categorical_accuracy(h, y))

        tfu.summary.scalar("cost", cross_entropy)
        tfu.summary.scalar("accuracy", accuracy)

        return dict(
            cost=cross_entropy,
            accuracy=accuracy,
        )

    file_writer = tf.summary.FileWriter(trial.file_path("summary"))

    train_summary = tfu.SummaryAccumulator()
    train_summary.add_file_writer(file_writer)

    valid_summary = tfu.SummaryAccumulator()
    valid_summary.add_file_writer(file_writer)

    updates = tfu.UpdatesAccumulator()

    with tfu.variable_scope("train"), train_summary, updates:
        train_out = model()
        tfu.updates.adam(train_out["cost"])

    with tfu.variable_scope("valid", deterministic=True), valid_summary:
        valid_out = model()

    # enable XLA
    config = tf.ConfigProto()
    config.graph_options.optimizer_options.global_jit_level = tf.OptimizerOptions.ON_1
    sess = tf.InteractiveSession(config=config)

    tfu.counter.set_session(sess)
    sess.run(tf.global_variables_initializer())

    # load previous iteration
    if 0:
        prev_trial = du.trial.TrialState(trial_name=trial.trial_name,
                                         iteration_num=9)
        tfu.serialization.load_variables(
            prev_trial.file_path("final_variables"))

    summary_printer = tfu.SummaryPrinter()
    summary_printer.add_recipe("trial_info", trial)
    summary_printer.add_recipe("progress")
    summary_printer.add_recipe("iter")
    summary_printer.add_recipe("time")
    summary_printer.add_recipe("s_per_iter")
    summary_printer.add_recipe("x min+iter",
                               "valid_epoch/cost",
                               format="%.4g")
    summary_printer.add_recipe("x max+iter",
                               "valid_epoch/accuracy",
                               format="%.4g")
    summary_printer.add_recipe("add_finals",
                               ["train_epoch/cost",
                                "train_epoch/accuracy",
                                "valid_epoch/cost",
                                "valid_epoch/accuracy"],
                               format="%.4g")
    train_summary.add_summary_printer(summary_printer)
    valid_summary.add_summary_printer(summary_printer)

    train_fn = tfu.tf_fn(sess=sess,
                         inputs={"x": x,
                                 "y": y},
                         outputs=train_out,
                         ops=[updates, train_summary, tfu.counter.step_op()])
    train_fn = tfu.wrap.output_nan_guard(train_fn)
    train_fn = tfu.wrap.split_input(train_fn, split_size=BATCH_SIZE)
    train_fn = tfu.wrap.shuffle_input(train_fn)
    train_fn = tfu.wrap.format_output_keys(train_fn, "train_epoch/%s")
    train_fn = tfu.wrap.update_summary_printer(train_fn, summary_printer)

    valid_fn = tfu.tf_fn(sess=sess,
                         inputs={"x": x,
                                 "y": y},
                         outputs=valid_out,
                         ops=[valid_summary])
    valid_fn = tfu.wrap.split_input(valid_fn, split_size=BATCH_SIZE)
    valid_fn = tfu.wrap.format_output_keys(valid_fn, "valid_epoch/%s")
    valid_fn = tfu.wrap.update_summary_printer(valid_fn, summary_printer)

    while tfu.counter.get_count_value() < expected_count:
        with du.timer("epoch"):
            train_res = train_fn(train)
            valid_res = valid_fn(valid)
            print(summary_printer.to_org_list())

    if 0:
        tfu.serialization.dump_variables(trial.file_path("final_variables"))


text:
import numpy as np
import scipy

import sys

from .. import utilits as ut

from ._simsig_tools import _check_list,_rand_uniform

__all__ = ['generator_base']

from abc import ABC, abstractmethod
class generator_base(ABC):
    _signal_parameters = [] # dictionary list
    _system_parameters = {}
    _signal_parameters_dict_default = {}
    _system_parameters_dcit_default = {}

    #------------------------------------------------------------------------------------
    @abstractmethod
    def __init__(self, *args, **kwargs):
        self.set_signal_parameters_dict_default()

    #------------------------------------------------------------------------------------
    @abstractmethod
    def set_system_parameters(self, *args, **kwargs):
        pass

    #------------------------------------------------------------------------------------
    #@prepared method
    def get_system_parameters(self):
        '''
        Get systemparameters.
        '''
        return self._system_parameters

    #------------------------------------------------------------------------------------
    #@prepared method
    def set_signal_parameters_dict_default(self):
        '''
        Set dictionary list by default dictionary.
        '''
        self.clear_signal_parameters()
        self.add_signal_parameters_dict_default()

    #------------------------------------------------------------------------------------
    #@prepared method
    def clear_signal_parameters(self):
        '''
        Clear dictionary list.
        '''
        self._signal_parameters = []

    #------------------------------------------------------------------------------------
    #@prepared method
    def add_signal_parameters_dict_default(self):
        '''
        Add default dictionary to the dictionary list.
        '''
        self._signal_parameters.append(self._signal_parameters_dict_default.copy())

    #------------------------------------------------------------------------------------
    #@prepared method
    def add_signal_parameters_dicts(self, dictionaries = None ):
        '''
        Add dictionaries to the dictionary list.

        Parameters
        -----------
        * dictionaries: dict,
             dictionaries for list.

        Notes
        -----------
        * list, np.ndarray and tuple formats of input are applied.

        '''
        if(dictionaries):
            if(type(dictionaries) in [list, tuple, np.ndarray] ):

                #TODO: add check check format
                for dict_i in dictionaries:
                    self._signal_parameters.append( dict_i )
        else:
            self.add_signal_parameters_dict_default()

    #------------------------------------------------------------------------------------
    #@prepared method
    def set_signal_parameters_dicts(self, dictionaries = None ):
        '''
        Set dictionary list.
        Notes: list, np.ndarray and tuple formats of input are applied.
        param dictionaries: - dictionaries for list.
        '''
        self.clear_signal_parameters()
        self.add_signal_parameters_dicts(dictionaries)

    #------------------------------------------------------------------------------------
    #@prepared method
    def remove_signal_parameters(self, n_of_parameters_dict = -1):
        '''
        Remove one dictionary form dictionary list.

        Parameters
        ------------
        n_of_parameters_dict: dict,
            number of the dict to remove.
        '''
        if(n_of_parameters_dict>=len(self._signal_parameters)):raise(IndexError)
        else:self._signal_parameters.pop(n_of_parameters_dict)

    #------------------------------------------------------------------------------------
    @abstractmethod
    def make_signal_parameters_dict(self, *args, **kwargs):
        '''
        Make signal_parameters dictionary.
        '''
        signal_parameters_dict = self.get_signal_parameters_dict_default()
        return signal_parameters_dict

    #------------------------------------------------------------------------------------
    @abstractmethod
    def add_signal_parameters(self, *args, **kwargs):
        '''
        Add signal_parameters.
        Notes
        --------
        * formats of the input: float, list, tuple, ndarray.
        '''
        pass
    #------------------------------------------------------------------------------------
    @abstractmethod
    def set_signal_parameters(self, *args, **kwargs):
        '''
        Set signal_parameters.
        Notes
        --------
        * formats of the input: float, list, tuple, ndarray.
        '''
        pass

    #------------------------------------------------------------------------------------
    @abstractmethod
    def add_random_signal_parameters(self, *args, **kwargs):
        '''
        Add random uniformly distributed signal_parameters.

        Notes
        -----------
        * callbacks doesnot applied for this function.
        '''
        pass

    #------------------------------------------------------------------------------------
    #@prepared method
    def get_signal_parameters_dicts(self):
        '''
        Get signal parameters dictionary list.
        '''
        return self._signal_parameters

    #------------------------------------------------------------------------------------
    #@prepared method
    def print_signal_parameters(self ):
        '''
        Print signal parameters dictionary list.
        '''
        print(30*'-')
        for i,dict_i in enumerate(self._signal_parameters): print(i,dict_i)

    #------------------------------------------------------------------------------------
    #@prepared method
    def get_signal_parameters_dict_default(self ):
        '''
        Get signal parameters default dictionary.
        '''
        return self._signal_parameters_dict_default.copy()

    #------------------------------------------------------------------------------------
    #@prepared method
    def _sim_signal(self, snr_db=None, signal_indexes = None):
        '''
        Simulate signal in form corresponing the class.

        Parameters
        ------------
        * snr_db: float,
            is the signal-to-noise ratio (in dB) for whole signal.
        * signal_indexes: int,
            index if signals in the signal dictionary.

        Notes
        ---------
        * int,float,complex, list, np.ndarray and tuple
           formats of signal_indexes are applied.
        * SNR in dB.
        * If callback is not None:
          s = callback(s) (format of callback = f(x)),
          if callback is None it does not applied.

        '''
        params = np.array(self._signal_parameters)

        if signal_indexes is None:
            signal_indexes = np.arange(params.shape[0])

        elif (isinstance(signal_indexes,(int,float,complex))):
            signal_indexes = np.append([],signal_indexes)

        elif (type(signal_indexes) not in [list, tuple,np.ndarray]):
            raise ValueError(' applied types of signal_indexes are: int,')

        if (np.max(signal_indexes)>params.shape[0]):
            raise ValueError(' indexs range exeed length of parameters, ',params.shape[0])

        params = params[signal_indexes]

        for i,param_dict in enumerate(params):

            signew = self._sim_one_sig(param_dict)
            if(i==0):
                self._signal = signew
            else:
                self._signal = np.add(self._signal,signew )

        if(snr_db):
            self._signal = ut.awgn(self._signal, snr_db, units = 'db')

    #------------------------------------------------------------------------------------
    #@prepared method
    def get_signal(self, snr_db=None):
        '''
        Get signal in form corresponing to the class.

        Parameters
        ------------
        * snr_db: float,
            is the signal-to-noise ratio
            (in dB) for whole signal.

        Returns
        -----------
        * signal: 1d ndarray.

        Notes
        -----------

        '''
        self._sim_signal(snr_db=snr_db, signal_indexes = None)
        return self._signal

    #------------------------------------------------------------------------------------
    def __call__(self, snr_db=None):
        '''
        Get signal in form corresponing to the class.

        Parameters
        ------------
        * snr_db: float,
            is the signal-to-noise ratio
            (in dB) for whole signal.

        Returns
        -----------
        * signal: 1d ndarray.

        Notes
        -----------

        '''
        return self.get_signal(snr_db)

    #------------------------------------------------------------------------------------
    @abstractmethod
    def _sim_one_sig(self, sig_param):
        '''
        Simulate harmonic (actually fmcw signal).

        Parameters
        -----------
        * sig_param: dict,
            dictionary of signal parameters.

        Returns
        ----------
        * signal:1d naarray,
            simulated signal.
        '''
        sig = np.zeros([],dtype = np.complex)

        if (callback in  ['None', None]):
              return sig

        elif type(callback ) is not list:
            callback = list([callback])

            for callback_i in callback:
                sig = callback_i(sig)

            return sig

ext: py

max_stars_repo_name: diarmaidocualain/dsatools

max_issues_repo_name: Jimmy-INL/dsatools

max_forks_repo_name: Jimmy-INL/dsatools

reason: import numpy,import scipy


text:
"""
costs according to supply systems
"""

import numpy as np
import pandas as pd
from geopandas import GeoDataFrame as gpdf
import itertools
import cea.config
import cea.inputlocator
from cea.analysis.costs.equations import calc_capex_annualized, calc_opex_annualized

__author__ = "Jimeno A. Fonseca"
__copyright__ = "Copyright 2020, Architecture and Building Systems - ETH Zurich"
__credits__ = ["Jimeno A. Fonseca", " Emanuel Riegelbauer"]
__license__ = "MIT"
__version__ = "0.1"
__maintainer__ = "Daren Thomas"
__email__ = "cea@arch.ethz.ch"
__status__ = "Production"

def costs_main(locator, config):
    # get local variables
    capital = config.costs.capital
    operational = config.costs.operational

    # get demand
    demand = pd.read_csv(locator.get_total_demand())

    # get the databases for each main system
    cooling_db, hot_water_db, electricity_db, heating_db = get_databases(demand, locator)

    # COSTS DUE TO HEATING SERIVICES (EXCEPT HOTWATER)
    heating_final_services = ['OIL_hs', 'NG_hs', 'WOOD_hs', 'COAL_hs', 'GRID_hs', 'DH_hs']
    costs_heating_services_dict = calc_costs_per_energy_service(heating_db, heating_final_services)

    # COSTS DUE TO HOT WATER SERVICES
    hot_water_final_services = ['OIL_ww', 'NG_ww', 'WOOD_ww', 'COAL_ww', 'GRID_ww', 'DH_ww']
    costs_hot_water_services_dict = calc_costs_per_energy_service(hot_water_db, hot_water_final_services)

    # COSTS DUE TO COOLING SERVICES
    cooling_final_services = ['GRID_cs', 'GRID_cdata', 'GRID_cre', 'DC_cs']
    costs_cooling_services_dict = calc_costs_per_energy_service(cooling_db, cooling_final_services)

    # COSTS DUE TO ELECTRICITY SERVICES
    electricity_final_services = ['GRID_pro', 'GRID_l', 'GRID_aux', 'GRID_v', 'GRID_a', 'GRID_data', 'GRID_ve']
    costs_electricity_services_dict = calc_costs_per_energy_service(electricity_db, electricity_final_services)

    # COMBINE INTO ONE DICT
    result = dict(itertools.chain(costs_heating_services_dict.items(), costs_hot_water_services_dict.items(),
                                  costs_cooling_services_dict.items(), costs_electricity_services_dict.items()))

    # sum up for all fields
    # create a dict to map from the convention of fields to the final variables
    mapping_dict = {'_capex_total_USD': 'Capex_total_sys_USD',
                    '_opex_fixed_USD': 'Opex_fixed_sys_USD',
                    '_opex_var_USD': 'Opex_var_sys_USD',
                    '_opex_USD': 'Opex_sys_USD',
                    # all system annualized
                    '_capex_a_USD': 'Capex_a_sys_USD',
                    '_opex_a_var_USD': 'Opex_a_var_sys_USD',
                    '_opex_a_fixed_USD': 'Opex_a_fixed_sys_USD',
                    '_opex_a_USD': 'Opex_a_sys_USD',
                    '_TAC_USD': 'TAC_sys_USD',
                    # building_scale_systems
                    '_capex_total_building_scale_USD': 'Capex_total_sys_building_scale_USD',
                    '_opex_building_scale_USD': 'Opex_sys_building_scale_USD',
                    '_capex_a_building_scale_USD': 'Capex_a_sys_building_scale_USD',
                    '_opex_a_building_scale_USD': 'Opex_a_sys_building_scale_USD',
                    # district_scale_systems
                    '_capex_total_district_scale_USD': 'Capex_total_sys_district_scale_USD',
                    '_opex_district_scale_USD': 'Opex_sys_district_scale_USD',
                    '_capex_a_district_scale_USD': 'Capex_a_sys_district_scale_USD',
                    '_opex_a_district_scale_USD': 'Opex_a_sys_district_scale_USD',
                    # city_scale_systems
                    '_capex_total_city_scale_USD': 'Capex_total_sys_city_scale_USD',
                    '_opex_city_scale_USD': 'Opex_sys_city_scale_USD',
                    '_capex_a_city_scale_USD': 'Capex_a_sys_city_scale_USD',
                    '_opex_a_city_scale_USD': 'Opex_a_sys_city_scale_USD',
                    }
    # initialize the names of the variables in the result to zero
    n_buildings = demand.shape[0]
    for _, value in mapping_dict.items():
        result[value] = np.zeros(n_buildings)

    # loop inside the results and sum the results
    for field in result.keys():
        for key, value in mapping_dict.items():
            if key in field:
                result[value] += result[field]

    # add name and create dataframe
    result.update({'Name': demand.Name.values})
    result_out = pd.DataFrame(result)

    # save dataframe
    result_out.to_csv(locator.get_costs_operation_file(), index=False, float_format='%.2f',  na_rep='nan')

def calc_costs_per_energy_service(database, heating_services):
    result = {}
    for service in heating_services:
        # TOTALS
        result[service + '_capex_total_USD'] = (database[service + '0_kW'].values *
                                                database['efficiency'].values *  # because it is based on the end use
                                                database['CAPEX_USD2015kW'].values)

        result[service + '_opex_fixed_USD'] = (result[service + '_capex_total_USD'] * database['O&M_%'].values / 100)

        result[service + '_opex_var_USD'] = database[service + '_MWhyr'].values * database[
            'Opex_var_buy_USD2015kWh'].values * 1000

        result[service + '_opex_USD'] = result[service + '_opex_fixed_USD'] + result[service + '_opex_var_USD']

        # ANNUALIZED
        result[service + '_capex_a_USD'] = np.vectorize(calc_capex_annualized)(result[service + '_capex_total_USD'],
                                                                               database['IR_%'],
                                                                               database['LT_yr'])

        result[service + '_opex_a_fixed_USD'] = np.vectorize(calc_opex_annualized)(result[service + '_opex_fixed_USD'],
                                                                                database['IR_%'],
                                                                                database['LT_yr'])

        result[service + '_opex_a_var_USD'] = np.vectorize(calc_opex_annualized)(result[service + '_opex_var_USD'],
                                                                                database['IR_%'],
                                                                                database['LT_yr'])

        result[service + '_opex_a_USD'] = np.vectorize(calc_opex_annualized)(result[service + '_opex_USD'],
                                                                             database['IR_%'],
                                                                             database['LT_yr'])

        result[service + '_TAC_USD'] = result[service + '_opex_a_USD'] + result[service + '_capex_a_USD']

        # GET CONNECTED AND DISCONNECTED
        for field in ['_capex_total_USD', '_capex_a_USD', '_opex_USD', '_opex_a_USD']:
            field_district = field.split("_USD")[0] + "_district_scale_USD"
            field_building_scale = field.split("_USD")[0] + "_building_scale_USD"
            field_city_scale = field.split("_USD")[0] + "_city_scale_USD"
            result[service + field_district], \
            result[service + field_building_scale], \
            result[service + field_city_scale] = np.vectorize(calc_scale_costs)(result[service + field],
                                                                                database['scale'])
    return result

def calc_scale_costs(value, flag_scale):
    if flag_scale == "BUILDING":
        district = 0.0
        building = value
        city = 0.0
    elif flag_scale == "DISTRICT":
        district = value
        building = 0.0
        city = 0.0
    elif flag_scale == "CITY":
        district = 0.0
        building = 0.0
        city = value
    elif flag_scale == "NONE":
        if value == 0.0 or np.isnan(value):
            district = 0.0
            building = 0.0
            city = 0.0
        else:
            raise ValueError("the scale is NONE but somehow there is a cost here?"
                             " the inputs of SUPPLY database may be wrong")
    else:
        raise ValueError("the scale in the system is {}, this is not a valid argument"
                         "valid arguments are CITY, DISTRICT, BUILDING, NONE".format(flag_scale))
    return district, building, city

def get_databases(demand, locator):
    supply_systems = gpdf.from_file(locator.get_building_supply()).drop('geometry', axis=1)
    data_all_in_one_systems = pd.read_excel(locator.get_database_supply_assemblies(), sheet_name=None)
    factors_heating = data_all_in_one_systems['HEATING']
    factors_dhw = data_all_in_one_systems['HOT_WATER']
    factors_cooling = data_all_in_one_systems['COOLING']
    factors_electricity = data_all_in_one_systems['ELECTRICITY']
    factors_resources = pd.read_excel(locator.get_database_feedstocks(), sheet_name=None)
    # get the mean of all values for this
    factors_resources_simple = [(name, values['Opex_var_buy_USD2015kWh'].mean()) for name, values in
                                factors_resources.items()]
    factors_resources_simple = pd.DataFrame(factors_resources_simple,
                                            columns=['code', 'Opex_var_buy_USD2015kWh']).append(
        # append NONE choice with zero values
        {'code': 'NONE'}, ignore_index=True).fillna(0)
    # local variables
    # calculate the total operational non-renewable primary energy demand and CO2 emissions
    ## create data frame for each type of end use energy containing the type of supply system use, the final energy
    ## demand and the primary energy and emissions factors for each corresponding type of supply system
    heating_costs = factors_heating.merge(factors_resources_simple, left_on='feedstock', right_on='code')[
        ['code_x', 'feedstock', 'scale', 'efficiency', 'Opex_var_buy_USD2015kWh', 'CAPEX_USD2015kW', 'LT_yr', 'O&M_%',
         'IR_%']]
    cooling_costs = factors_cooling.merge(factors_resources_simple, left_on='feedstock', right_on='code')[
        ['code_x', 'feedstock', 'scale', 'efficiency', 'Opex_var_buy_USD2015kWh', 'CAPEX_USD2015kW', 'LT_yr', 'O&M_%',
         'IR_%']]
    dhw_costs = factors_dhw.merge(factors_resources_simple, left_on='feedstock', right_on='code')[
        ['code_x', 'feedstock', 'scale', 'efficiency', 'Opex_var_buy_USD2015kWh', 'CAPEX_USD2015kW', 'LT_yr', 'O&M_%',
         'IR_%']]
    electricity_costs = factors_electricity.merge(factors_resources_simple, left_on='feedstock', right_on='code')[
        ['code_x', 'feedstock', 'scale', 'efficiency', 'Opex_var_buy_USD2015kWh', 'CAPEX_USD2015kW', 'LT_yr', 'O&M_%',
         'IR_%']]
    heating = supply_systems.merge(demand, on='Name').merge(heating_costs, left_on='type_hs', right_on='code_x')
    dhw = supply_systems.merge(demand, on='Name').merge(dhw_costs, left_on='type_dhw', right_on='code_x')
    cooling = supply_systems.merge(demand, on='Name').merge(cooling_costs, left_on='type_cs', right_on='code_x')
    electricity = supply_systems.merge(demand, on='Name').merge(electricity_costs, left_on='type_el', right_on='code_x')
    return cooling, dhw, electricity, heating

def main(config):
    locator = cea.inputlocator.InputLocator(scenario=config.scenario)

    print('Running system-costs with scenario = %s' % config.scenario)

    costs_main(locator=locator, config=config)

if __name__ == '__main__':
    main(cea.config.Configuration())

reason: import numpy


text:
"""
Fit data for transit alone, after "detrending" the stellar variability.
"""
import os
import numpy as np, pandas as pd, matplotlib.pyplot as plt, pymc3 as pm
from os.path import join

from timmy.modelfitter import ModelFitter, ModelParser
import timmy.plotting as tp
from timmy.convenience import get_clean_rv_data
from timmy.priors import initialize_prior_d

from timmy.paths import RESULTSDIR

def main(modelid):

    assert modelid == 'rv'

    make_threadsafe = 0

    phaseplot = 0
    cornerplot = 1

    OVERWRITE = 1
    REALID = 'TOI_837'

    PLOTDIR = os.path.join(
        RESULTSDIR, '{}_{}_fitting_results'.format(REALID, modelid)
    )
    if not os.path.exists(PLOTDIR):
        os.mkdir(PLOTDIR)
    PLOTDIR = os.path.join(PLOTDIR, '20200525')
    if not os.path.exists(PLOTDIR):
        os.mkdir(PLOTDIR)

    print(42*'#')
    print(modelid)
    print(42*'#')

    pklpath = os.path.join(
        os.path.expanduser('~'), 'local', 'timmy',
        '{}_model_{}.pkl'.format(REALID, modelid)
    )
    np.random.seed(42)

    rv_df = get_clean_rv_data()

    mp = ModelParser(modelid)

    prior_d = initialize_prior_d(mp.modelcomponents)

    m = ModelFitter(
        modelid, rv_df, prior_d, plotdir=PLOTDIR, pklpath=pklpath,
        overwrite=OVERWRITE
    )

    var_names = ['r_star', 'logg_star', 'log_K', 'period', 'ecc', 'omega', 'phase',
                 't0', 'S_tot', 'ell', 'means', 'sigmas']

    print(pm.summary(m.trace, var_names=var_names))
    summdf = pm.summary(m.trace, var_names=var_names, round_to=10,
                        kind='stats', stat_funcs={'median':np.nanmedian},
                        extend=True)

    if make_threadsafe:
        pass

    else:
        if cornerplot:
            outpath = join(PLOTDIR, '{}_{}_cornerplot.png'.format(REALID, modelid))
            tp.plot_cornerplot(prior_d, m, outpath)

        import IPython; IPython.embed()
        #FIXME: save the SAMPLES from the trace (rather than the trace itself,
        #b/c of the local functions definined within the model context)

        assert 0
        #FIXME plots below

        if phaseplot:
            outpath = join(PLOTDIR, '{}_{}_phaseplot.png'.format(REALID, modelid))
            tp.plot_phasefold(m, summdf, outpath)

        assert 0

        if sampleplot:
            outpath = join(PLOTDIR, '{}_{}_sampleplot.png'.format(REALID, modelid))
            tp.plot_sampleplot(m, outpath, N_samples=100)

if __name__ == "__main__":
    main('rv')

max_stars_repo_name: lgbouma/timmy

max_issues_repo_name: lgbouma/timmy

max_forks_repo_name: lgbouma/timmy


text:
"""
Preprocess a raw json dataset into hdf5/json files for use in data_loader.lua

Input: json file that has the form
[{ file_path: 'path/img.jpg', captions: ['a caption', ...] }, ...]
example element in this list would look like
{'captions': [u'A man with a red helmet on a small moped on a dirt road. ', u'Man riding a motor bike on a dirt road on the countryside.', u'A man riding on the back of a motorcycle.', u'A dirt path with a young person on a motor bike rests to the foreground of a verdant area with a bridge and a background of cloud-wreathed mountains. ', u'A man in a red shirt and a red hat is on a motorcycle on a hill side.'], 'file_path': u'val2014/COCO_val2014_000000391895.jpg', 'id': 391895}

This script reads this json, does some basic preprocessing on the captions
(e.g. lowercase, etc.), creates a special UNK token, and encodes everything to arrays

Output: a json file and an hdf5 file
The hdf5 file contains several fields:
/images is (N,3,256,256) uint8 array of raw image data in RGB format
/labels is (M,max_length) uint32 array of encoded labels, zero padded
/label_start_ix and /label_end_ix are (N,) uint32 arrays of pointers to the
  first and last indices (in range 1..M) of labels for each image
/label_length stores the length of the sequence for each of the M sequences

The json file has a dict that contains:
- an 'ix_to_word' field storing the vocab in form {ix:'word'}, where ix is 1-indexed
- an 'images' field that is a list holding auxiliary information for each image,
  such as in particular the 'split' it was assigned to.
"""

### source activate py2.7 ###

from __future__ import absolute_import
from __future__ import division
from __future__ import print_function

import os
import json
import sys
import re
import argparse
from random import shuffle, seed
import string
# non-standard dependencies:
import h5py
import numpy as np
import skimage.io
from PIL import Image

sys.path.append('/home/yz9244/Up-Down-Captioner/external/coco/PythonAPI/')
from pycocotools.coco import COCO

VizWiz_ANN_PATH = '/home/yz9244/Up-Down-Captioner/bottom-up-attention/data/VizWiz/annotations/'
COCO_TRAIN_VOCAB_PATH = '/home/yz9244/AoANet/data/cocotalk_vocab.json'

corrupt_list = [37093]

SENTENCE_SPLIT_REGEX = re.compile(r'(\W+)') # Split on any non-alphanumeric character

def split_sentence(sentence):
  """ break sentence into a list of words and punctuation """
  toks = []
  for word in [s.strip().lower() for s in SENTENCE_SPLIT_REGEX.split(sentence.strip()) if len(s.strip()) > 0]:
    # Break up any words containing punctuation only, e.g. '!?', unless it is multiple full stops e.g. '..'
    if all(c in string.punctuation for c in word) and not all(c in '.' for c in word):
      toks += list(word)
    else:
      toks.append(word)
  # Remove '.' from the end of the sentence -
  # this is EOS token that will be populated by data layer
  if toks[-1] != '.':
    return toks
  return toks[:-1]

def load_caption_vocab(vocab_path=COCO_TRAIN_VOCAB_PATH):
  info = json.load(open(vocab_path))
  ix_to_word = info['ix_to_word']
  vocab = []
  vocab_size = len(ix_to_word)
  for i in range(vocab_size):
    vocab.append(ix_to_word[str(i+1)])
  return vocab

def build_vocab(params, data_split, base_vocab):
  vocabulary = load_caption_vocab(vocab_path=base_vocab)
  offset = len(vocabulary)
  print("number of words in the base vocab: %d\n" % (offset))
  return vocabulary

def encode_captions(params, data_split, wtoi):
  """
  encode all captions into one large array, which will be 1-indexed.
  also produces label_start_ix and label_end_ix which store 1-indexed
  and inclusive (Lua-style) pointers to the first and last caption for
  each image in the dataset.
  """
  max_length = params['max_length']
  img_count = 0
  caption_count = 0
  for dataset in data_split:
    annFile='%s/VizWiz_Captions_v1_%s.json' % (VizWiz_ANN_PATH, dataset)
    coco = COCO(annFile)
    for image_id,anns in coco.imgToAnns.iteritems():
      if image_id in corrupt_list:
        continue
      img_count += 1
      caption_count += len(anns)
  N = img_count
  M = caption_count # total number of captions

  label_arrays = []
  label_start_ix = np.zeros(N, dtype='uint32') # note: these will be one-indexed
  label_end_ix = np.zeros(N, dtype='uint32')
  label_length = np.zeros(M, dtype='uint32')
  caption_counter = 0
  counter = 1
  img_counter = 0
  imgInfo = []

  for dataset in data_split:
    annFile='%s/VizWiz_Captions_v1_%s.json' % (VizWiz_ANN_PATH, dataset)
    coco = COCO(annFile)
    for image_id,anns in coco.imgToAnns.iteritems():
      if image_id in corrupt_list:
        continue
      image_info = coco.imgs[image_id]
      #image_path = '%s/%s' % (image_info['file_name'].split('_')[1], image_info['file_name'])
      image_path = '%s' % (image_info['file_name'])
      jimg = {}
      if dataset == 'train' or dataset == 'val':
        jimg['split'] = 'train'
      elif dataset == 'test':
        jimg['split'] = 'test'
      jimg['file_path'] = image_path
      jimg['id'] = image_info['id']
      n = len(anns)
      assert n > 0, 'error: some image has no captions'
      Li = np.zeros((n, max_length), dtype='uint32')
      for j, ann in enumerate(anns):
        caption_sequence = split_sentence(ann['caption'])
        label_length[caption_counter] = min(max_length, len(caption_sequence))
        caption_counter += 1
        for k,w in enumerate(caption_sequence):
          if k < max_length:
            if not w in wtoi:
              w = 'UNK'
            Li[j,k] = wtoi[w]
      # note: word indices are 1-indexed, and captions are padded with zeros
      label_arrays.append(Li)
      label_start_ix[img_counter] = counter
      label_end_ix[img_counter] = counter + n - 1
      img_counter += 1
      counter += n
      imgInfo.append(jimg)

  L = np.concatenate(label_arrays, axis=0) # put all the labels together
  assert L.shape[0] == M, 'lengths don\'t match? that\'s weird'
  assert np.all(label_length > 0), 'error: some caption had no words?'

  print('encoded captions to array of size ', L.shape)
  return N, L, label_start_ix, label_end_ix, label_length, imgInfo

def main(params):
  seed(123) # make reproducible

  # create the vocab integrating MSCOCO
  vocab = build_vocab(params, ['train', 'val'], base_vocab=COCO_TRAIN_VOCAB_PATH)
  itow = {i+1:w for i,w in enumerate(vocab)} # a 1-indexed vocab translation table
  wtoi = {w:i+1 for i,w in enumerate(vocab)} # inverse table
  # encode captions in large arrays, ready to ship to hdf5 file
  N, L, label_start_ix, label_end_ix, label_length, imgInfo = encode_captions(params, ['train', 'val', 'test'], wtoi)
  # create output h5 file
  f_lb = h5py.File(params['output_h5']+'_pretrained_label.h5', "w")
  f_lb.create_dataset("labels", dtype='uint32', data=L)
  f_lb.create_dataset("label_start_ix", dtype='uint32', data=label_start_ix)
  f_lb.create_dataset("label_end_ix", dtype='uint32', data=label_end_ix)
  f_lb.create_dataset("label_length", dtype='uint32', data=label_length)
  f_lb.close()
  # create output json file
  out = {}
  out['ix_to_word'] = itow # encode the (1-indexed) vocab
  out['images'] = imgInfo
  json.dump(out, open(params['output_json']+'_pretrained.json', 'w'))
  print('wrote ', params['output_json']+'_pretrained.json')

if __name__ == "__main__":

  parser = argparse.ArgumentParser()

  # input json
  parser.add_argument('--output_json', default='data/vizwiztalk', help='output json file')
  parser.add_argument('--output_h5', default='data/vizwiztalk', help='output h5 file')
  parser.add_argument('--images_root', default='', help='root location in which images are stored, to be prepended to file_path in input json')

  # options
  parser.add_argument('--max_length', default=16, type=int, help='max length of a caption, in number of words. captions longer than this get clipped.')
  parser.add_argument('--word_count_threshold', default=5, type=int, help='only words that occur more than this number of times will be put in vocab')

  args = parser.parse_args()
  params = vars(args) # convert to ordinary dict
  print('parsed input parameters:')
  print(json.dumps(params, indent = 2))
  main(params)

max_stars_repo_name: Yinan-Zhao/AoANet

max_issues_repo_name: Yinan-Zhao/AoANet

max_forks_repo_name: Yinan-Zhao/AoANet


text:
import numpy as np
import pandas as pd
import math
import json
from matplotlib import pyplot as plt
import seaborn as sns

from datetime import datetime

import matplotlib.pyplot as plt
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter
from mpl_toolkits.mplot3d import Axes3D
from kmodes.kprototypes import KPrototypes

from scipy import interpolate

import sqlite3
from sqlalchemy import create_engine

engine = create_engine('sqlite:///data/data.db')

# read in the json files
portfolio = pd.read_json('data/portfolio.json', orient='records', lines=True)
profile = pd.read_json('data/profile.json', orient='records', lines=True)
transcript = pd.read_json('data/transcript.json', orient='records', lines=True)

transcript.rename (columns = {'time': 'time in hours'}, inplace = True)

def interquartile_range (data, cols, x):

    '''
        In Data Science Interquartile range is defined as the area shown below
        upper = 75 % + (75 % - 25 %)
        lower = 25 % - (75 % - 25 %)

        values outside this range are propably false and do influence the analyzes in a wrong direction
    '''

    for column in cols:
        upper = data.describe () [column]['75%'] + ((data.describe () [column]['75%'] - data.describe () [column]['25%'])*x)
        lower = data.describe () [column]['25%'] - ((data.describe () [column]['75%'] - data.describe () [column]['25%'])*x)
        data = data [data [column] < upper]
        data = data [data [column] > lower]
    return data

def ETL_portfolio ():

    '''
        ETL Part of the Portfolio DF:

        - Duration is in days, matching the transcripts hours
        - the channel col is kind of a dictionary, so I get the dummie Variables.

    '''

    portfolio ['durationHOURS']= portfolio['duration']*24

    i = 0

    channel_dummies = []
    channel_dummies = pd.DataFrame(channel_dummies)
    while i < len (portfolio):
        channel_dummies = channel_dummies.append (pd.DataFrame(pd.get_dummies(portfolio ['channels'].loc [i]).sum (axis = 0)).T, sort = False)
        i += 1
    channel_dummies = channel_dummies.fillna (0).astype (int).reset_index ().drop (['index'], axis = 1)
    channel_dummies

    portfolio_a = portfolio.merge (channel_dummies, left_index = True, right_index = True)
    portfolio_a

    return portfolio_a

def get_Transcript_values (start, end):
    '''
        the transcript value col is also a dictionary
        for faster access it has to be rewritten

    '''

    value = []
    amount = []

    i = start
    while i < end:
        try:
            len (transcript['value'].loc[i].get ('offer id'))> 1
            if len (transcript['value'].loc[i].get ('offer id'))> 1:
                value.append (transcript['value'].loc[i].get ('offer id'))
                amount.append (0)

            elif len (transcript['value'].loc[i].get ('offer_id'))> 1:
                value.append (transcript['value'].loc[i].get ('offer_id'))
                amount.append (0)

            elif len (transcript['value'].loc[i].get ('amount'))> 1:
                value.append (0)
                amount.append (transcript['value'].loc[i].get ('amount'))
            else:
                value.append (0)
                amount.append (0)
        except:
            if (transcript['value'].loc[i].get ('amount')) is not None:
                value.append (0)
                amount.append (transcript['value'].loc[i].get ('amount'))

            elif len (transcript['value'].loc[i].get ('offer_id'))> 1:
                value.append (transcript['value'].loc[i].get ('offer_id'))
                amount.append (0)

            else:
                value.append (99999)
                amount.append (99999)

        i += 1

    value = pd.DataFrame (value)
    values = pd.get_dummies (value[0]).merge (value, left_index = True, right_index = True)
    amount = pd.DataFrame (amount)
    transcript_a = transcript.merge (amount, left_index = True, right_index = True)
    transcript_a = transcript_a.merge (values, left_index =True, right_index = True)
    transcript_a.rename (columns = {'0_y': 'offer_id', 0: 'amount'}, inplace = True)

    return transcript_a

def ETL_transcript (transcript):

    '''
        I do find it helpfull to also have the expiratioin time
        and the offer type within that data frame.
    '''

    expiration = []
    offer_type = []
    i = 0

    while i < len(transcript):

        if transcript.iloc [i] ['event'] == 'offer received':
            expiration.append (transcript.iloc [i] ['time in hours'] + int(portfolio [portfolio ['id'] == transcript.iloc [i] ['offer_id']]['durationHOURS']))
            offer_type.append (portfolio.set_index (['id']).loc [transcript.iloc [i] ['offer_id']] ['offer_type'])
        else:
            expiration.append (transcript.iloc [i] ['time in hours'])
            offer_type.append ('0')

        i += 1

    transcript ['expiration'] = expiration
    transcript ['offer_type'] = offer_type

    return transcript

def get_visit_data (transcript):

    '''
        Not quite sure if i can get a detailled user profil.
        Human behavior can be shown between weeks.
        I am therefor setting up a day dummie variable and a week variable

        be aware, that it is unclear wether day 0 of the week is a monday or maybe a thursday.
        But all days labeled day-0 are the same day in the week.
    '''

    day_of_month = (transcript['time in hours']/24).astype (int)
    day_of_month = pd.DataFrame (day_of_month)

    day_of_week = day_of_month - 7*(day_of_month / 7).astype (int)
    day_of_week = pd.DataFrame (day_of_week['time in hours'])

    week_of_month = (day_of_month / 7).astype (int)

    day_of_month.rename (columns = {'time in hours': 'DayOfMonth'}, inplace = True)
    day_of_week.rename (columns = {'time in hours' : 'DayOfWeek'}, inplace = True)
    week_of_month.rename (columns = {'time in hours' : 'WeekOfMonth'}, inplace = True)

    transcript = transcript.merge (day_of_week, left_index = True, right_index = True)
    transcript = transcript.merge (week_of_month, left_index = True, right_index = True)
    transcript.rename (columns = {'0_x': 'count', '0_y': 'offer_id'}, inplace = True)
    transcript = transcript.merge (day_of_month, left_index = True, right_index = True)

    return transcript

def split_on_event (data):

    '''
        Spliting the data set for better performance

    '''

    received = data [data ['event'] == 'offer received']

    viewed = data [data ['event'] == 'offer viewed']

    completed = data [data ['event'] == 'offer completed']

    transactions = data [data ['event'] == 'transaction']

    return received, viewed, completed, transactions

def ETL_profiles (profile):
    '''
        drop Nan and none gender, aswell as age == 118

    '''
    profile = profile.dropna (axis = 0, subset =  ['income'])
    genders = pd.get_dummies (profile ['gender'])
    profile = profile [profile ['age'] < 100]
    profile = profile.dropna (axis = 0, subset = ['income'])

    '''
        Member since is a date, but I rewrite it in weeks till the very moment.
    '''

    profiles_b = profile

    i = 0
    member_since = []

    while i < len (profiles_b):
        member_since.append (int(int(str(datetime.today () - datetime (int(profiles_b.iloc [i]['became_member_on'].astype (str)[0:4]),
                              int (profiles_b.iloc [i]['became_member_on'].astype (str)[4:6]),
                              int (profiles_b.iloc [i]['became_member_on'].astype (str)[6:]))).split (' ')[0])/ 7))
        i += 1

    genders = genders.merge (pd.DataFrame (member_since), left_index = True, right_index = True)

    profile = profile.merge (genders, left_index = True, right_index = True)
    profile.rename (columns = {0: 'member_x_weeks'}, inplace = True)

    return (profile)

'''
________________________________________________________________________________

'''

portfolio = ETL_portfolio()

print ('ETL Portfolio')

transcript = get_Transcript_values (0, len (transcript))
transcript = interquartile_range (transcript, ['amount'], 3)
transcript = ETL_transcript (transcript)
transcript = get_visit_data(transcript)

print ('ETL Transcript')

received, viewed, completed, transactions = split_on_event (transcript)

profile = ETL_profiles (profile)
print ('ETL Profiles')

'''
________________________________________________________________________________
'''

def get_user_interactions (user_id):
    user_int = transcript [transcript ['person'] == user_id]
    return user_int

def get_offer_response (i):

    # data = 3 Parts... + transactions
    received_index = i
    received_matrix = received.loc [i]
    offer_id = (received.loc [i] ['offer_id'])
    person = (received.loc [i] ['person'])
    received_time = received.loc [i] ['time in hours']
    expiration_time = received.loc [i] ['expiration']

    viewed_yn = []
    viewtoexpire = []
    completed_yn = []
    completed_yn_unknown = []
    compacttime = []
    comp_matrix = []
    liste = []

    view_matrix = viewed [viewed ['person'] == person][viewed [viewed ['person'] == person]['offer_id'] == offer_id]
    view_matrix = view_matrix [view_matrix ['time in hours'] >= received_time]
    view_matrix = view_matrix [view_matrix ['time in hours'] <= expiration_time]

    try:
        viewtime = view_matrix.iloc [0]['time in hours']
    except:
        viewtime = -98

    comp_matrix = completed [completed ['person'] == person][completed [completed ['person'] == person]['offer_id']== offer_id]
    comp_matrix = comp_matrix [comp_matrix ['time in hours'] >= received_time]

    try:
        completed_time = (comp_matrix.iloc [0]['time in hours'])
    except:
        completed_time = -99

    if viewtime >= 0:
        viewed_yn = 1
        viewtoexpire = expiration_time - viewtime
    else:
        viewed_yn = 0
        viewtoexpire = 0

    if completed_time >= 0:
        if viewtime >= 0:
            compacttime =  completed_time - viewtime
            completed_yn = 1
            completed_yn_unknown = 0
        else:
            compacttime = 0
            completed_yn_unknown =  1
            completed_yn = 0
    else:
        completed_yn_unknown = 0
        compacttime = 0
        completed_yn = 0

    index_received = received_index
    liste = index_received

    index_viewed = -1
    try:
        index_viewed = view_matrix.index [0].astype (int)
        liste.append (index_viewed)
    except: index_viewed = index_viewed

    index_completed = -1
    try:
        '''
            transaction Data, that leads into a completed offer is one record above!!
        '''
        index_completed = (comp_matrix.index [0].astype (int) - 1)

        if (transcript.loc [index_completed]['event'] == 'transaction'):
            index_completed = index_completed

        elif (transcript.loc [index_completed -1 ]['event'] == 'transaction'):
            index_completed = index_completed -1

        else:
            index_completed = -1

        liste.append (index_completed)
    except: index_completed = index_completed

    start = viewtime
    end = min (completed_time, expiration_time)

    visits = transactions [transactions ['person'] == person][transactions [transactions ['person'] == person] ['time in hours'] >= start]
    visits = visits [visits ['time in hours'] < end]
    visits_count = len (visits [visits ['time in hours'] < end])

    vector = ['offer_id', 'person', 'viewed_yn', 'viewtoexpire', 'completed_yn', 'compacttime', 'completed_yn_unknown', 'index_received', 'index_viewed', 'index_completed', 'visits_count']

    offer_response_matrix = pd.DataFrame ([offer_id, person, viewed_yn, viewtoexpire, completed_yn, compacttime, completed_yn_unknown, index_received, index_viewed, index_completed, visits_count])

    return offer_response_matrix

def get_user_transaction_data ():

    '''
        kind of a customer analyses, but keep in mind, that there is no split between control and test group.
        additional, since there is no control period, this data is distracted
    '''

    transcript_test = transcript

    visits_person = transactions.groupby (['person'])['amount'].value_counts().unstack().sum (axis = 1)

    frequent_visits = []
    frequent_visits = pd.DataFrame (frequent_visits)

    total_spent_each_person = pd.DataFrame(transcript_test [transcript_test ['event'] == 'transaction'].groupby (['person']).sum ())
    frequent_visits = frequent_visits.merge (total_spent_each_person, left_index = True, right_index = True, how = 'outer')

    value_person = transcript_test.groupby (['person'])['amount'].sum ().sort_values(ascending = False)
    deep_value_person = transcript_test [transcript_test ['event'] == 'transaction'].groupby (['person'])['amount'].describe ()
    i = 0
    std_arr = []

    while i < len (deep_value_person):
        try: std_arr.append (int(deep_value_person['std'].iloc[i]*100)/100)
        except: std_arr.append (-1)
        i+= 1

    deep_value_person = deep_value_person.reset_index ().merge (pd.DataFrame (std_arr), left_index = True, right_index = True).set_index (['person'])
    deep_value_person = deep_value_person.drop(['std'], axis = 1).rename (columns = {0: 'std'})

    profile_return = profile.merge (deep_value_person, left_on = 'id', right_on = 'person')

    return profile_return

# get the matrix (viewtime, ...) on each offer!

def get_reaction_matrix (list_of_customers):

    a = 0
    # count variable

    matrix = []
    matrix = pd.DataFrame (matrix)

    while a <  len (list_of_customers):
        '''
        If you want a smaller set of customers
        '''

        person = list_of_customers [a]
        offers = get_user_interactions (person) [get_user_interactions (person) ['event'] == 'offer received']
        i = 0
        while i < len (offers):
            n = offers.index [i]
            matrix = matrix.append (get_offer_response (n).T)
            i += 1
        a += 1

    vector = ['offer_id', 'person', 'viewed_yn', 'viewtoexpire', 'completed_yn', 'compacttime', 'completed_yn_unknown', 'index_received', 'index_viewed', 'index_completed', 'visits_count']

    matrix.columns = pd.DataFrame(vector).T.iloc [0]

    received_reaction_matrix = received.merge (matrix, left_index = True, right_on = 'index_received')

    matrix.columns = pd.DataFrame(vector).T.iloc [0]

    received_reaction_matrix = received.merge (matrix, left_index = True, right_on = 'index_received')
    return received_reaction_matrix

def store_df (data, name):

    '''
        store data function in SQL
    '''

    name_in = name

    data.to_sql (name_in, engine, if_exists = 'replace')
    print (name_in, ' stored')
    return

'''
________________________________________________________________________________
'''

profile = get_user_transaction_data ()
print ('ETL offer response Matrix')

list_of_customers = (list(profile ['id']))

received_reaction_matrix = get_reaction_matrix (list_of_customers)
received_reaction_matrix

visits_person = transactions.groupby (['person'])['amount'].value_counts().unstack().sum (axis = 1)

profiles_in_use = interquartile_range (profile, ['age', 'income', 'max'], 1)
profiles_in_use
print ('Profiles are cleaned')

'''
    time to store all data

'''

### portfolio:
name = 'portfolio'
store_df (portfolio.drop (['channels'], axis = 1), name)

### profile:
name = 'profiles'
store_df (profiles_in_use, name)

### transcript:
name = 'transcript'
store_df (transcript.drop (['value'], axis = 1).astype (str), name)

### received-reaction-matrix:
name = 'received_reaction_matrix'
store_df (received_reaction_matrix.drop (['value'], axis = 1).astype (str), name)

reason: import numpy,from scipy


text:
# coding: utf-8
#import networkx as nx
#import matplotlib.pyplot as plt
#import operator
#from collections import defaultdict
#from collections import Counter
#from collections import deque
# from functools import reduce
# from math import log
import copy
import operator
import re
import sys
import time
from itertools import product
from operator import add, mul

import numpy as np

operators = {
    '+': add,
    '*': mul,
}

def find_parens(s):
    ret = {}
    pstack = []

    for i, c in enumerate(s):
        if c == '(':
            pstack.append(i)
        elif c == ')':
            ret[pstack.pop()] = i

    return ret

def calculate(s, DBG=True):
    if DBG:
        print(s)

    parens = find_parens(s)

    lhs = 0
    ret = 0
    operator = '+'
    i = 0
    while i < len(s):
        if s[i].isdigit():
            lhs = int(s[i])
            ret = operators[operator](ret, lhs)
            i = i+1
        elif s[i] in operators.keys():
            operator = s[i]
            i = i+1
        elif s[i] == '(':
            ss = s[i+1:parens[i]]
            if DBG:
                print(ss)
            lhs = calculate(ss, DBG)
            ret = operators[operator](ret, lhs)
            i = parens[i] + 1
        elif s[i] == ' ':
            i = i+1
            pass
    return ret

def boom(input_val, DBG=True):

    ret = 0
    for expr in input_val:
        value = calculate(expr, DBG)
        ret = ret + value

    return ret

def test(cc=None, expected=None, DBG=False):
    start_millis = int(round(time.time() * 1000))
    result = boom(cc, DBG)
    stop_millis = int(round(time.time() * 1000))
    result = str(result)
    expected = str(expected)
    flag = (result == expected)
    if(expected == "None"):
        print("*** "+str(cc) + " *** -> Result = "+str(result))
    else:
        print("*** "+str(cc) + " *** -> Result = "+str(result) +
              " -> success = " + str(flag) + " -> expected " + expected)
    print((stop_millis-start_millis), "ms", int((stop_millis-start_millis) /
                                                1000), "s", int((stop_millis-start_millis)/1000/60), "min")
    return flag

t1 = """1 + 2 * 3 + 4 * 5 + 6"""
tt1 = t1.splitlines()
test(tt1, 71, True)
# sys.exit()

t1 = """1 + (2 * 3) + (4 * (5 + 6))"""
tt1 = t1.splitlines()
test(tt1, 51, True)
# sys.exit()

t1 = """2 * 3 + (4 * 5)"""
tt1 = t1.splitlines()
test(tt1, 26, True)
# sys.exit()
t1 = """5 + (8 * 3 + 9 + 3 * 4 * 3)"""
tt1 = t1.splitlines()
test(tt1, 437, True)
# sys.exit()
t1 = """5 * 9 * (7 * 3 * 3 + 9 * 3 + (8 + 6 * 4))"""
tt1 = t1.splitlines()
test(tt1, 12240, True)
# sys.exit()
t1 = """((2 + 4 * 9) * (6 + 9 * 8 + 6) + 6) + 2 + 4 * 2"""
tt1 = t1.splitlines()
test(tt1, 13632, True)
# sys.exit()

INPUT_FILE = "input-d18.txt"
f = open(INPUT_FILE, "r")
contents = f.read()
puzzle_input = contents.splitlines()
f.close()

ret = boom(puzzle_input, DBG=False)
print(ret)

# part 1 = 650217205854

max_stars_repo_name: jbudynek/advent-of-code

max_issues_repo_name: jbudynek/advent-of-code

max_forks_repo_name: jbudynek/advent-of-code

reason: import numpy,import networkx


text:
import pytest

import os
import numpy as np
import scipy.io

from pylgr import legendre_gauss_radau as LGR

TOL = 1e-10

# Reference data computed with Matlab scripts
test_data_path = os.path.join('unit_tests', 'test_data', 'LGR_diff.mat')
LGR_diff_mat = scipy.io.loadmat(test_data_path)
LGR_diff_reference = {}

# Process data to make it easier to deal with
for case_num in LGR_diff_mat:
    if case_num[0] == 'n':
        case_key = int(case_num[1:])
        LGR_diff_reference[case_key] = {}
        for key in ['tau','w','D']:
            LGR_diff_reference[case_key][key] = np.squeeze(
                LGR_diff_mat[case_num][key][0][0]
            )

@pytest.mark.parametrize('n', [-1,0,1,2])
@pytest.mark.parametrize('fun_name', ['make_LGR_nodes', 'make_LGR'])
def test_make_LGR_small_n(n, fun_name):
    fun = getattr(LGR, fun_name)
    with pytest.raises(ValueError):
        fun(n)

@pytest.mark.parametrize('n', [0,1,2])
@pytest.mark.parametrize(
    'fun_name', ['make_LGR_weights', 'make_LGR_diff_matrix']
)
def test_make_LGR_small_tau(n, fun_name):
    fun = getattr(LGR, fun_name)
    tau = np.random.rand(n)
    with pytest.raises(ValueError):
        fun(tau)

@pytest.mark.parametrize('n', LGR_diff_reference.keys())
def test_make_LGR(n):
    tau, w, D = LGR.make_LGR(n)
    assert np.allclose(tau, LGR_diff_reference[n]['tau'])
    assert np.allclose(w, LGR_diff_reference[n]['w'])
    assert np.allclose(D, LGR_diff_reference[n]['D'])

@pytest.mark.parametrize('n', np.arange(3,19,3))
def test_LGR_basic_int_diff(n):
    '''
    Tests some basic identities: w*x=0, Dx=ones, Dx**2=2x, w*(Dx)=2.
    '''
    tau, w, D = LGR.make_LGR(n)
    assert np.abs(np.dot(w, tau)) < TOL
    assert np.allclose(np.matmul(D, tau), np.ones_like(tau))
    assert np.allclose(np.matmul(D, tau**2), 2.*tau)
    assert np.isclose(np.dot(w, np.matmul(D, tau)), 2.)

@pytest.mark.parametrize('n', np.arange(3,19,3))
def test_LGR_integrate(n):
    '''
    LGR should be able to integrate a polynomial of degree 2n-2 to machine
    precision.
    '''
    # Generate a random polynomial of degree 2n-2
    degree = 2*n - 2
    coef = np.random.randn(degree+1)
    P = np.polynomial.polynomial.Polynomial(coef)
    expected_integral = P.integ(lbnd=-1.)(1.)

    tau = LGR.make_LGR_nodes(n)
    w = LGR.make_LGR_weights(tau)
    LGR_integral = np.dot(w, P(tau))

    assert np.isclose(LGR_integral, expected_integral)

@pytest.mark.parametrize('n', np.arange(3,19,3))
def test_LGR_differentiate(n):
    '''
    LGR should be able to differentiate a polynomial of degree n-1 to machine
    precision.
    '''
    # Generate a random polynomial of degree n-1
    degree = n - 1
    coef = np.random.randn(degree+1)
    P = np.polynomial.polynomial.Polynomial(coef)

    tau, w, D = LGR.make_LGR(n)
    LGR_derivative = np.matmul(D, P(tau))

    assert np.allclose(LGR_derivative, P.deriv()(tau))

@pytest.mark.parametrize('n_dims', [1,2,3])
def test_LGR_multivariate_integrate(n_dims):
    '''
    Test integration of degree 2n-2 polymomials in n_dims dimensions.
    '''
    n = 10

    # Generate random polynomial of degree 2n-2
    degree = 2*n - 2
    coef = np.random.randn(n_dims, degree+1)
    P = [np.polynomial.polynomial.Polynomial(coef[d]) for d in range(n_dims)]
    expected_integral = [P[d].integ(lbnd=-1.)(1.) for d in range(n_dims)]

    tau = LGR.make_LGR_nodes(n)
    w = LGR.make_LGR_weights(tau)
    P_mat = np.vstack([P[d](tau) for d in range(n_dims)])
    LGR_integral = np.matmul(P_mat, w)

    assert np.allclose(LGR_integral, expected_integral)

@pytest.mark.parametrize('n_dims', [1,2,3])
def test_LGR_multivariate_differentiate(n_dims):
    '''
    Test differentiation of degree n-1 polymomials in n_dims dimensions.
    '''
    n = 10

    # Generate a random polynomial of degree n-1
    degree = n - 1
    coef = np.random.randn(n_dims, degree+1)
    P = [np.polynomial.polynomial.Polynomial(coef[d]) for d in range(n_dims)]

    tau, w, D = LGR.make_LGR(n)
    P_mat = np.vstack([P[d](tau) for d in range(n_dims)])

    expected_derivative = [P[d].deriv()(tau) for d in range(n_dims)]
    LGR_derivative = np.matmul(P_mat, D.T)

    assert np.allclose(LGR_derivative, expected_derivative)

max_stars_repo_name: Tenavi/PyLGR

max_issues_repo_name: Tenavi/PyLGR

max_forks_repo_name: Tenavi/PyLGR


text:
import numpy as np
import pandas as pd
from sklearn.decomposition import PCA, FastICA, KernelPCA, TruncatedSVD, IncrementalPCA, LatentDirichletAllocation, MiniBatchDictionaryLearning, MiniBatchSparsePCA, NMF, SparsePCA, SparseCoder, FactorAnalysis

class factorAnalysis(FactorAnalysis):
    feature_name = 'factoranalysis'
    def transform(self, X):
        data = super().transform(X)
        cols = [f'{self.feature_name}_{i}' for i in range(data.shape[1])]
        return pd.DataFrame(data, columns = cols, index = X.index)

class pCA(PCA):
    feature_name = 'pca'
    def transform(self, X):
        data = super().transform(X)
        cols = [f'{self.feature_name}_{i}' for i in range(data.shape[1])]
        return pd.DataFrame(data, columns = cols, index = X.index)

class fastICA(FastICA):
    feature_name = 'fastica'
    def transform(self, X):
        data = super().transform(X)
        cols = [f'{self.feature_name}_{i}' for i in range(data.shape[1])]
        return pd.DataFrame(data, columns = cols, index = X.index)

class incrementalPCA(IncrementalPCA):
    feature_name = 'incrementalpca '
    def transform(self, X):
        data = super().transform(X)
        cols = [f'{self.feature_name}_{i}' for i in range(data.shape[1])]
        return pd.DataFrame(data, columns = cols, index = X.index)

class kernelPCA(KernelPCA):
    feature_name = 'kernelpca'
    def transform(self, X):
        data = super().transform(X)
        cols = [f'{self.feature_name}_{i}' for i in range(data.shape[1])]
        return pd.DataFrame(data, columns = cols, index = X.index)

class latentDirichletAllocation(LatentDirichletAllocation):
    feature_name = 'incrementalpca'
    def transform(self, X):
        data = super().transform(X)
        cols = [f'{self.feature_name}_{i}' for i in range(data.shape[1])]
        return pd.DataFrame(data, columns = cols, index = X.index)

class miniBatchDictionaryLearning(MiniBatchDictionaryLearning):
    feature_name = 'minibatchdictionarylearning'
    def transform(self, X):
        data = super().transform(X)
        cols = [f'{self.feature_name}_{i}' for i in range(data.shape[1])]
        return pd.DataFrame(data, columns = cols, index = X.index)

class miniBatchSparsePCA(MiniBatchSparsePCA):
    feature_name = 'minibatchsparsepca'
    def transform(self, X):
        data = super().transform(X)
        cols = [f'{self.feature_name}_{i}' for i in range(data.shape[1])]
        return pd.DataFrame(data, columns = cols, index = X.index)

class nMF(NMF):
    feature_name = 'nmf'
    def transform(self, X):
        data = super().transform(X)
        cols = [f'{self.feature_name}_{i}' for i in range(data.shape[1])]
        return pd.DataFrame(data, columns = cols, index = X.index)

class sparsePCA(SparsePCA):
    feature_name = 'sparsepca'
    def transform(self, X):
        data = super().transform(X)
        cols = [f'{self.feature_name}_{i}' for i in range(data.shape[1])]
        return pd.DataFrame(data, columns = cols, index = X.index)

class sparseCoder(SparseCoder):
    feature_name = 'sparsecoder'
    def transform(self, X):
        data = super().transform(X)
        cols = [f'{self.feature_name}_{i}' for i in range(data.shape[1])]
        return pd.DataFrame(data, columns = cols, index = X.index)

class truncatedSVD(TruncatedSVD):
    feature_name = 'truncatedsvd'
    def transform(self, X):
        data = super().transform(X)
        cols = [f'{self.feature_name}_{i}' for i in range(data.shape[1])]
        return pd.DataFrame(data, columns = cols, index = X.index)


text:
#!/usr/bin/env python3

import pandas

# import matplotlib.pyplot as plt

import numpy as np
from scipy.optimize import curve_fit, minimize

def sigmoid(value):
  return 1.0 / (1.0 + np.exp(-value))

def max_point(xdata, ydata):
  variables = len(xdata)
  weights = (variables + 1) * (variables + 2) // 2

  def fit_func(X, *params):
    prod = 1.0
    assert len(X) == variables
    ret = params[0]
    var_id = 1
    for i in range(variables):
      x = X[i]
      ret += x * params[var_id]
      var_id += 1
      for j in range(i, variables):
        y = X[j]
        ret += x * y * params[var_id]
        var_id += 1
    assert var_id == len(params)
    ret = sigmoid(ret / 10.0)
    print(ret)
    return ret

  p0 = np.array((1,) * weights)
  popt, pcov = curve_fit(fit_func, xdata, ydata, p0)
  print('popt: {}'.format(popt))

  minmax = (xdata.min(1), xdata.max(1))
  bounds = np.array(list(zip(*minmax)))
  print('minmax:\n', minmax)
  print('bounds:\n', bounds)
  fm = lambda x: -fit_func(x, *popt)
  r = minimize(fm, bounds.mean(1), bounds=bounds)
  return r["x"], fit_func(r["x"], *popt)

def main():

  data = pandas.read_csv('train_result.csv')

  print(data.describe())

  buckets = {}
  count = {}

  color_map = {
    -1: (1, 0, 0),
     0: (0, 0, 1),
     1: (0, 1, 1),
  }

  column_names = [
          'up_weight',
          'down_weight',
          'left_weight',
          'right_weight',
  ]

  # X = []
  # Y = []
  # c = []
  # c = np.array([color_map[v] for v in data.result])

  print('Samples: {}'.format(len(data)))

  result = np.array((data.result + 1) / 2)

  columns = np.array([data[name] for name in column_names]) / 10000

  maximum_coords, max_value = max_point(columns, result)

  print('Maximum value: {}'.format(max_value))
  for name, value in zip(column_names, maximum_coords):
    print('{}: {}'.format(name, value))

  # buckets = {}
  # bucket_count = {}
  # for i in range(len(data)):
  #   k = data.approach_weapon[i]
  #   v = result[i]
  #   if not k in buckets:
  #     buckets[k] = 0
  #     bucket_count[k] = 0
  #   buckets[k] += v
  #   bucket_count[k] += 1

  # for k in buckets:
  #   buckets[k] /= bucket_count[k]

  # print(sorted(buckets.items()))

if __name__ == '__main__':
  main()

max_stars_repo_name: rahenri/sawbot

max_issues_repo_name: rahenri/sawbot

max_forks_repo_name: rahenri/sawbot


text:
import sys,os
import numpy as np
import pandas as pd
import seaborn as sns
sns.set(style='whitegrid')
import matplotlib.pyplot as plt

def main(argv):
  outfile = argv[0]
  metric = ['Precision','Precision','Precision','Recall','Recall','Recall','F-score','F-score','F-score']
  method = ['Heuristic','Random Forests','Deep learning','Heuristic','Random Forests','Deep learning','Heuristic','Random Forests','Deep learning']
  values = [69.23, 71.10, 72.18, 68.68, 74.53, 64.41, 68.94, 72.57, 66.96]
  df = pd.DataFrame({'metric': metric, 'method': method, 'values': values})
  ax = sns.barplot(x='metric', y='values', hue='method', data=df)
  ax.set(xlabel='Performance metric', ylabel='Values (%)')
  plt.legend(loc=2, bbox_to_anchor=(1.05, 1), borderaxespad=0.)
  plt.title('Wake vs. Sleep classification on AMC dataset (105 subjects)')
  plt.ylim([0,100])

  def plot_val(metric, st_idx):
    ax.text(metric-0.4,values[st_idx]+2,str(values[st_idx]))
    ax.text(metric-0.1,values[st_idx+1]+2,str(values[st_idx+1]))
    ax.text(metric+0.2,values[st_idx+2]+2,str(values[st_idx+2]))
  plot_val(0,0)
  plot_val(1,3)
  plot_val(2,6)

  fig = ax.get_figure()
  fig.savefig(outfile, bbox_inches='tight')

if __name__ == '__main__':
  main(sys.argv[1:])

max_stars_repo_name: wadpac/SleepStageClassification

max_issues_repo_name: wadpac/SleepStageClassification

max_forks_repo_name: wadpac/SleepStageClassification


text:
import numpy as np
from numpy.random import randn
import unittest
from niwqg import CoupledModel
from niwqg import QGModel

class QGNIWTester(unittest.TestCase):
    """ A class for testing the QGNIW kernel (real and complex 2d ffts)
        Note: 1d plane wave pass test with machine precision
              2d (slanted) plane wave has an error O(10^{-13})
                                                                        """
    def setUp(self):
        self.m =  CoupledModel.Model(use_filter=False,nu4=1e14,nu4w=0.)
        self.m.tmax = 10*self.m.dt
        k, l  = 2*np.pi*5/self.m.L, 2*np.pi*9/self.m.L
        self.qi = np.sin(k*self.m.x + l*self.m.y)
        self.m.set_q(self.qi)
        self.m.set_phi(self.qi*0)

    def test_hyperviscosity(self, rtol=1e-15):
        """ Test if the hyperviscosity implementation simply damps
                the Fourier coefficiants individualy. """

        self.m.run()
        qfh = self.m.fft(self.qi)*np.exp(-self.m.nu4*self.m.wv4*self.m.tmax)
        self.assertTrue(np.allclose(qfh,self.m.qh,rtol=rtol), 'Implementation of\
                             hypterviscosity is broken in CoupledModel')

class QGTester(unittest.TestCase):
    """ A class for testing the QG model (rffts)
            Note: 1d plane wave pass test with machine precision
                  2d (slanted) plane wave has an error O(10^{-13})
                                                                        """
    def setUp(self):
        self.m =  QGModel.Model(use_filter=False, nu4 = 1e14)
        self.m.tmax = 100*self.m.dt
        k, l  = 2*np.pi*5/self.m.L, 2*np.pi*9/self.m.L
        self.qi = np.sin(k*self.m.x + l*self.m.x)
        self.m.set_q(self.qi)

    def test_hyperviscosity(self, rtol=1e-15):
        """ Test if the hyperviscosity implementation simply damps
                the Fourier coefficiants individualy. """

        self.m.run()
        qfh = self.m.fft(self.qi)*np.exp(-self.m.nu4*self.m.wv4*self.m.tmax)
        self.assertTrue(np.allclose(qfh,self.m.qh,rtol=rtol), 'Implementation of\
                             hypterviscosity is broken QGModel')

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

reason: import numpy,from numpy


text:
import numpy as np

class EleMomentumRemi():
    def __init__(self,t0_ns=0,x0_mm=0,y0_mm=0,vjetx_mmPns=0,vjety_mmPns=0,
    l_mm=None,U_V=None):
        self.t0 = t0_ns
        self.x0 = x0_mm
        self.y0 = y0_mm
        self.vjetx = vjetx_mmPns
        self.vjety = vjety_mmPns
        self.U = U_V
        self.l = l_mm

        self.mmPns2au = 0.4571028904957055
        self.VPmm2mmPns = 1.75882014896e-1
        self.amu2au = 1836.15
        self.au2tesla = 2.35e5
        self.au2mm = 5.28e-8

    def CalcPtr(self,B_tesla,omega_mns,x_mm,y_mm,t_ns):
        R = np.sqrt((x_mm-self.x0)**2 + (y_mm-self.y0)**2)
        angle = omega_mns*(t_ns-self.t0)
        angle %= 2*np.pi
        Ptr = (B_tesla/self.au2tesla)*R/self.au2mm/(2*np.abs(np.sin(angle/2)))

        theta = (np.arctan2(y_mm-self.y0,x_mm-self.x0)+2*np.pi)%(2*np.pi)
        phi = theta - angle/2
        Px = Ptr*np.cos(phi) - self.mmPns2au*self.vjetx
        Py = Ptr*np.sin(phi) - self.mmPns2au*self.vjety

        return Ptr, phi, Px, Py

    def CalcR(self,B_tesla,omega_mns,Ptr,t_ns):
        angle = omega_mns*(t_ns-self.t0)
        angle %= 2*np.pi
        R = self.au2mm*(2*np.abs(np.sin(angle/2)))*Ptr/(B_tesla/self.au2tesla)
        return R

    def CalcPzOneAcc(self,m_amu,t_ns):
        return self.amu2au*m_amu*self.mmPns2au*self.l/(t_ns-self.t0) - 8.04e-2*self.U*(t_ns-self.t0)/(2*self.l)

    def CalcPzOneAccApprox(self,ta_ns,ta0_ns,sfc=None):
        if sfc is not None:
            return sfc*(ta0_ns-ta_ns)
        else:
            return 8.04e-2*self.U*(ta0_ns-ta_ns)/self.l


text:
from manimlib.imports import *
import numpy as np

class ODEField(GraphScene):
    CONFIG = {
        "x_min" : -4.5,
        "x_max" : 4.5,
        "y_min" : -3,
        "y_max" : 3,
        "graph_origin" : ORIGIN ,
        "function_color" : BLUE ,
        "axes_color" : GREEN,
        "center_point" : 0,
    }
    def construct(self):
        self.setup_axes(animate=False)
        def getField(xPos,yPos):
            magnitude = 0.5
            derivative = yPos-xPos
            angle = np.arctan(derivative)
            y = np.sin(angle)
            x = np.cos(angle)

            return magnitude*(RIGHT*x + UP*y)

        vectors = [Vector(getField(x,y), color = YELLOW).shift(RIGHT*x+UP*y)
        for x in np.arange(-10,10,1)
        for y in np.arange(-9,9,1)
        ] #List of vectors pointing to each grid point

        draw_field = VGroup(*vectors) #Pass list of vectors to create a VGroup

        self.play(ShowCreation(draw_field)) #Draw VGroup on screen
        self.wait(2)
        minX = -10
        maxX=10
        def solution(c):
            return lambda x: c*np.exp(x)+x+1

        def getUpdateFunction(startValue, endValue, f):

            def update(graph, alpha):
                currValue = interpolate(startValue, endValue, alpha)
                newGraph = self.get_graph(f(currValue),graph.color,x_max=maxX,x_min = minX)
                graph.become(newGraph)
                return graph

            return update
        fGraph = self.get_graph(solution(-1),BLUE,x_max=maxX,x_min = minX)
        self.play(ShowCreation(fGraph))
        self.play(UpdateFromAlphaFunc(fGraph,getUpdateFunction(-1,2,solution)),run_time=10)
        self.wait(3)
        dot = Dot(self.coords_to_point(0,1.25), color=YELLOW)
        self.play(ShowCreation(dot))
        self.play(UpdateFromAlphaFunc(fGraph,getUpdateFunction(2,0.25,solution)),run_time=3)
        #self.play(ShowCreation(exponential),run_time=4)
        self.wait(10)

class ConstantDerivativeV2(GraphScene):
    CONFIG = {
        "x_min" : -4.5,
        "x_max" : 4.5,
        "y_min" : -3,
        "y_max" : 3,
        "graph_origin" : ORIGIN ,
        "function_color" : BLUE ,
        "axes_color" : GREEN,
        "center_point" : 0,
    }
    def construct(self):
        self.setup_axes(animate=False)
        b=0.4
        def getField(xPos,yPos):
            magnitude = 0.5
            derivative = b
            angle = np.arctan(derivative)
            y = np.sin(angle)
            x = np.cos(angle)

            return magnitude*(RIGHT*x + UP*y)

        vectors = [Vector(getField(x,y), color = YELLOW).shift(RIGHT*x+UP*y)
        for x in np.arange(-10,10,1)
        for y in np.arange(-9,9,1)
        ] #List of vectors pointing to each grid point

        draw_field = VGroup(*vectors) #Pass list of vectors to create a VGroup

        self.play(ShowCreation(draw_field)) #Draw VGroup on screen
        self.wait(2)
        minX = -10
        maxX=10
        def solution(c):
            return lambda x: b*x+c
        def getLineGraph(c):
            return self.get_graph(solution(c),BLUE,x_max=maxX,x_min = minX)
        correctLine = getLineGraph(-1)
        lines = VGroup(*[getLineGraph(c) for c in np.arange(-5,5,0.5)])
        self.play(ShowCreation(correctLine), ShowCreation(lines))
        self.wait(3)
        dot = Dot(self.coords_to_point(0,-1), color=YELLOW)
        self.play(ShowCreation(dot))
        self.play(FadeOut(lines))
        #self.play(ShowCreation(exponential),run_time=4)
        self.wait(10)

class ConstantDerivative(GraphScene):
    CONFIG = {
        "x_min" : -4.5,
        "x_max" : 4.5,
        "y_min" : -3,
        "y_max" : 3,
        "graph_origin" : ORIGIN ,
        "function_color" : BLUE ,
        "axes_color" : GREEN,
        "center_point" : 0,
    }
    def construct(self):
        self.setup_axes(animate=False)
        b=0.4
        def getField(xPos,yPos):
            magnitude = 0.5
            derivative = b
            angle = np.arctan(derivative)
            y = np.sin(angle)
            x = np.cos(angle)

            return magnitude*(RIGHT*x + UP*y)

        vectors = [Vector(getField(x,y), color = YELLOW).shift(RIGHT*x+UP*y)
        for x in np.arange(-10,10,1)
        for y in np.arange(-9,9,1)
        ] #List of vectors pointing to each grid point

        draw_field = VGroup(*vectors) #Pass list of vectors to create a VGroup

        self.play(ShowCreation(draw_field)) #Draw VGroup on screen
        self.wait(2)
        minX = -10
        maxX=10
        def solution(c):
            return lambda x: b*x+c

        def getUpdateFunction(startValue, endValue, f):

            def update(graph, alpha):
                currValue = interpolate(startValue, endValue, alpha)
                newGraph = self.get_graph(f(currValue),graph.color,x_max=maxX,x_min = minX)
                graph.become(newGraph)
                return graph

            return update
        fGraph = self.get_graph(solution(-1),BLUE,x_max=maxX,x_min = minX)
        self.play(ShowCreation(fGraph))
        self.play(UpdateFromAlphaFunc(fGraph,getUpdateFunction(-1,2,solution)),run_time=10)
        self.wait(3)
        dot = Dot(self.coords_to_point(0,1.25), color=YELLOW)
        self.play(ShowCreation(dot))
        self.play(UpdateFromAlphaFunc(fGraph,getUpdateFunction(2,0.25,solution)),run_time=3)
        #self.play(ShowCreation(exponential),run_time=4)
        self.wait(10)

max_stars_repo_name: Prdcc/manim-Student-Shapers

max_issues_repo_name: Prdcc/manim-Student-Shapers

max_forks_repo_name: Prdcc/manim-Student-Shapers


text:
##############################################################################
#
# Copyright (c) 2009-2018 by The University of Queensland
# http://www.uq.edu.au
#
# Primary Business: Queensland, Australia
# Licensed under the Apache License, version 2.0
# http://www.apache.org/licenses/LICENSE-2.0
#
# Development until 2012 by Earth Systems Science Computational Center (ESSCC)
# Development 2012-2013 by School of Earth Sciences
# Development from 2014 by Centre for Geoscience Computing (GeoComp)
#
##############################################################################

from __future__ import division, print_function
__copyright__="""Copyright (c) 2009-2018 by The University of Queensland
http://www.uq.edu.au
Primary Business: Queensland, Australia"""
__license__="""Licensed under the Apache License, version 2.0
http://www.apache.org/licenses/LICENSE-2.0"""
__url__="https://launchpad.net/escript-finley"

"""
Author: Antony Hallam antony.hallam@uqconnect.edu.au
"""
############################################################FILE HEADER
# example06.py
# Create a 3 block fault and overburden style model in 2d using pycad
# meshing tools.

#######################################################EXTERNAL MODULES
import matplotlib
matplotlib.use('agg') #It's just here for automated testing
from esys.pycad import *
from esys.pycad.gmsh import Design
from esys.escript import *
import numpy as np
import pylab as pl #Plotting package
from cblib import toRegGrid, subsample, HAVE_NATGRID
from esys.escript.unitsSI import *
from esys.escript.linearPDEs import LinearPDE
import os, sys

try:
    # This imports the rectangle domain function
    from esys.finley import MakeDomain#Converter for escript
    HAVE_FINLEY = True
except ImportError:
    print("Finley module not available")
    HAVE_FINLEY = False
########################################################MPI WORLD CHECK
if getMPISizeWorld() > 1:
        import sys
        print("This example will not run in an MPI world.")
        sys.exit(0)

if not HAVE_NATGRID:
    print("This example requires that natgrid is available to matplotlib")

if HAVE_FINLEY and HAVE_NATGRID:
    #################################################ESTABLISHING VARIABLES
    # where to put output files
    save_path= os.path.join("data","example06")
    mkDir(save_path)

    Ttop=20*Celsius # temperature at the top
    qin=300.*Milli*W/(m*m) #our heat source temperature is now zero

    ################################################ESTABLISHING PARAMETERS
    # Overall Domain
    p0=Point(0.0,        0.0, 0.0)
    p1=Point(0.0,    -6000.0, 0.0)
    p2=Point(5000.0, -6000.0, 0.0)
    p3=Point(5000.0,     0.0, 0.0)

    ####################################################DOMAIN CONSTRUCTION
    l01=Line(p0, p1)
    l12=Line(p1, p2)
    l23=Line(p2, p3)
    l30=Line(p3, p0)

    # Generate Material Boundary
    p4=Point(0.0,    -2400.0, 0.0)
    p5=Point(2000.0, -2400.0, 0.0)
    p6=Point(3000.0, -6000.0, 0.0)
    p7=Point(5000.0, -2400.0, 0.0)

    # Create TOP BLOCK
    tbl1=Line(p0,p4)
    tbl2=Line(p4,p5)
    tbl3=Line(p5,p7)
    tbl4=Line(p7,p3)
    tblockloop = CurveLoop(tbl1,tbl2,tbl3,tbl4,l30)
    tblock = PlaneSurface(tblockloop)

    # Create BOTTOM BLOCK LEFT
    bbll1=Line(p4,p1)
    bbll2=Line(p1,p6)
    bbll3=Line(p6,p5)
    bbll4=-tbl2
    bblockloopl = CurveLoop(bbll1,bbll2,bbll3,bbll4)
    bblockl = PlaneSurface(bblockloopl)

    # Create BOTTOM BLOCK RIGHT
    bbrl1=Line(p6,p2)
    bbrl2=Line(p2,p7)
    bbrl3=-tbl3
    bbrl4=-bbll3
    bblockloopr = CurveLoop(bbrl1,bbrl2,bbrl3,bbrl4)
    bblockr = PlaneSurface(bblockloopr)

    #############################################EXPORTING MESH FOR ESCRIPT
    # Create a Design which can make the mesh
    d=Design(dim=2, element_size=200)
    # Add the subdomains and flux boundaries.
    d.addItems(PropertySet("top",tblock),\
                           PropertySet("bottomleft",bblockl),\
                           PropertySet("bottomright",bblockr),\
                           PropertySet("linebottom",bbll2, bbrl1))

    # Create the geometry, mesh and Escript domain
    d.setScriptFileName(os.path.join(save_path,"example06.geo"))

    d.setMeshFileName(os.path.join(save_path,"example06.msh"))
    domain=MakeDomain(d)
    print("Domain has been generated ...")

    # set up kappa (thermal conductivity across domain) using tags
    kappa=Scalar(0,Function(domain))
    kappa.setTaggedValue("top",2.0)
    kappa.setTaggedValue("bottomleft",10.0)
    kappa.setTaggedValue("bottomright",6.0)
    ##############################################################SOLVE PDE
    mypde=LinearPDE(domain)
    mypde.getSolverOptions().setVerbosityOn()
    mypde.setSymmetryOn()
    mypde.setValue(A=kappa*kronecker(domain))
    x=Solution(domain).getX()
    mypde.setValue(q=whereZero(x[1]-sup(x[1])),r=Ttop)
    qS=Scalar(0,FunctionOnBoundary(domain))
    qS.setTaggedValue("linebottom",qin)
    mypde.setValue(y=qS)
    print("PDE has been generated ...")
    ###########################################################GET SOLUTION
    T=mypde.getSolution()
    print("PDE has been solved ...")
    ###############################################################PLOTTING
    # show temperature:
    xi, yi, zi = toRegGrid(T, nx=50, ny=50)
    CS = pl.contour(xi,yi,zi,5,linewidths=0.5,colors='k')
    pl.clabel(CS, inline=1, fontsize=8)
    # show sub domains:
    tpg=np.array([p.getCoordinates() for p in tblockloop.getPolygon() ])
    pl.fill(tpg[:,0],tpg[:,1],'brown',label='2 W/m/k',zorder=-1000)
    bpgr=np.array([p.getCoordinates() for p in bblockloopr.getPolygon() ])
    pl.fill(bpgr[:,0],bpgr[:,1],'orange',label='6 W/m/k',zorder=-1000)
    bpgl=np.array([p.getCoordinates() for p in bblockloopl.getPolygon() ])
    pl.fill(bpgl[:,0],bpgl[:,1],'red',label='10 W/m/k',zorder=-1000)
    # show flux:
    xflux, flux=subsample(-kappa*grad(T), nx=20, ny=20)
    pl.quiver(xflux[:,0],xflux[:,1],flux[:,0],flux[:,1], angles='xy',color="white")
    # create plot
    pl.title("Heat Refraction across an anisotropic structure\n with flux vectors.")
    pl.xlabel("Horizontal Displacement (m)")
    pl.ylabel("Depth (m)")
    pl.legend()
    pl.savefig(os.path.join(save_path,"flux.png"))
    print("Flux has been plotted  ...")


text:
# -*- coding: utf-8 -*-
"""
Created on Fri Apr 10 09:43:32 2020

@author: muril
"""

def model_comp_plot(df_plot,
                    x_nm,
                    y_nm,
                    fv_nm,
                    hue_nm,
                    fn_out,
                    save_fig = True,
                    v_sub    = None,
                    p_sub    = ["#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7"],
                    size_p   = 4,
                    t_fsize  = 20,
                    t_ypos   = 0.93,
                    dpi_fig  = 300):

    #-------------------------------------------------------#
    #-------------------- model_comp_plot ------------------#
    #-------------------------------------------------------#
    #--- Goal:
    #---    Plot Simulated (lines) and Observed (markers) data against time
    #--- Parameters:
    #--- 	df_plot 	: dataframe with simulated and observed data for one or more variables
    #--- 	x_nm    	: Name of df_plot column to use as X axis
    #--- 	y_nm    	: Name of df_plot column to use as Y axis
    #--- 	fv_nm   	: Name of df_plot column to be split in different facets
    #---    hue_nm      : Name of df_plot column to be used as hue
    #--- 	fn_out  	: Filename of output figure
    #--- 	save_fig	: Save Figure (True/False)
    #--- 	size_p  	: Size of facet
    #--- 	t_fsize 	: Title fontsize
    #--- 	t_ypos  	: Title vertical position
    #--- 	dpi_fig 	: Figure Resolution
    #--- Author:
    #---    Murilo Vianna (murilodsv@gmail.com)
    #-------------------------------------------------------#

    import matplotlib.pyplot as plt
    import numpy as np
    from matplotlib.ticker import FormatStrFormatter

    #--- Get sim and obs data
    df_obs = df_plot[:][df_plot['type'] == 'obs']
    df_sim = df_plot[:][df_plot['type'] == 'sim']

    #--- get unique var and units
    df_uvar = df_plot[[fv_nm,'units']].drop_duplicates()
    v_nm = df_uvar[fv_nm].unique()

    #--- attribute unique colors to each hue
    hue_cols = {}
    c_i      = 0
    for h in df_plot[hue_nm].unique():
        hue_cols[h] = p_sub[c_i]
        c_i +=1
        if (c_i+1) > len(p_sub): c_i = 0

    #--- compute squared grid
    nv = len(v_nm)
    if np.sqrt(nv) - int(np.sqrt(nv)) == 0:
        nr_f = int(np.sqrt(nv))
        nc_f = int(np.sqrt(nv))
    else:
        nc_f = int(np.sqrt(nv))+1
        nr_f = int(nv / nc_f) + 1

    if nr_f * nc_f - nv == np.sqrt(nr_f * nc_f):
        nr_f = nr_f - 1

    x_size = nc_f   * size_p
    y_size = nr_f   * size_p

    #--- Initilize plot grid
    fig, axs = plt.subplots(nr_f, nc_f, figsize=(x_size, y_size))

    v = 0
    add_leg = True
    l_leg   = []
    for r_f in range(0, nr_f):
        for c_f in range(0,nc_f):

            #--- get axis obj
            if (nr_f * nc_f) == 1:      axs_rc = axs
            elif (nr_f * nc_f) == 2:    axs_rc = axs[c_f]
            else:                       axs_rc = axs[r_f, c_f]

            if (v+1) <= len(v_nm):

                #--- get variable
                var = v_nm[v]

                #--- get units
                u   = df_uvar['units'][df_uvar[fv_nm] == var].values[0]

                #--- plot simulations as lines
                for h in df_sim[hue_nm].unique():
                    f = (df_sim[fv_nm] == var) & (df_sim[hue_nm] == h)
                    axs_rc.plot(df_sim[x_nm][f],
                                df_sim[y_nm][f],
                                c = hue_cols[h],
                                linewidth=1)

                #--- plot observations as markers
                for h in df_obs[hue_nm].unique():

                    l = None
                    if add_leg:
                        l = h
                        l_leg = l_leg + [l]

                    f = (df_obs[fv_nm] == var) & (df_obs[hue_nm] == h)
                    axs_rc.scatter(df_obs[x_nm][f],
                                   df_obs[y_nm][f],
                                   edgecolors='black',
                                   facecolors=hue_cols[h],
                                   linewidths= 0.5,
                                   label = l)

                if add_leg and set(l_leg).issubset(df_obs[hue_nm]):
                    #axs[r_f, c_f].legend()
                    #fig.legend(loc = 'lower center', ncol = 4)
                    add_leg = False

                #--- Axis
                axs_rc.set(xlabel = x_nm,
                           ylabel = var+' ('+str(u)+')')

                #--- get limits
                max_vals = []
                min_vals = []
                if len(df_sim[y_nm][df_sim[fv_nm] == var].values) > 0:
                    max_vals = max_vals + [max(df_sim[y_nm][df_sim[fv_nm] == var].values)]
                    min_vals = min_vals + [min(df_sim[y_nm][df_sim[fv_nm] == var].values)]

                if len(df_obs[y_nm][df_obs[fv_nm] == var].values) > 0:
                    max_vals = max_vals + [max(df_obs[y_nm][df_obs[fv_nm] == var].values)]
                    min_vals = min_vals + [min(df_obs[y_nm][df_obs[fv_nm] == var].values)]

                max_y = max(max_vals)
                min_y = min(min_vals)

                if max_y - min_y < 0.01:
                    m_point = np.mean([max_y,min_y])
                    axs_rc.ylim([m_point-0.01, m_point+0.01])

                if max_y < 10:
                    axs_rc.yaxis.set_major_formatter(FormatStrFormatter('%.2f'))
                elif max_y < 100:
                    axs_rc.yaxis.set_major_formatter(FormatStrFormatter('%.1f'))
                elif max_y < 1000:
                    axs_rc.yaxis.set_major_formatter(FormatStrFormatter('%.0f'))

            else:
                #--- delete blank axis
                fig.delaxes(axs_rc)

            v += 1

    #fig.tight_layout()
    if (nr_f * nc_f) <= 2:
        fig.legend(loc = 'upper center', ncol = 4)
    else:
        fig.legend(loc = 'lower center', ncol = 4)

    #--- save results
    if save_fig:
        plt.savefig(fn_out+'.comp_plot.png', dpi = dpi_fig)

    #plt.show()
    plt.close()

max_stars_repo_name: Murilodsv/py-jules

max_issues_repo_name: Murilodsv/py-jules

max_forks_repo_name: Murilodsv/py-jules


text:
#!/usr/bin/env python
import zmq
import time
import dill
import logging
import threading
import numpy as np
import ami.comm
from ami import LogConfig
from ami.export.nt import NTBytes, NTObject, NTGraph, NTStore
from p4p.nt import NTScalar, NTNDArray
from p4p.server import Server, StaticProvider
from p4p.server.thread import SharedPV
from p4p.rpc import rpc, NTURIDispatcher
from p4p.util import ThreadedWorkQueue

logger = logging.getLogger(LogConfig.get_package_name(__name__))

class PvaExportPutHandler:

    def __init__(self, put=None, rpc=None):
        self._put = put
        self._rpc = rpc

    def put(self, pv, op):
        if self._put is not None:
            self._put(pv, op)

    def rpc(self, pv, op):
        if self._rpc is not None:
            self._rpc(pv, op)

class PvaExportRpcHandler:
    def __init__(self, ctx, addr):
        self.ctx = ctx
        self.addr = addr
        self.comms = {}

    def _get_comm(self, graph):
        if graph not in self.comms:
            self.comms[graph] = ami.comm.GraphCommHandler(graph, self.addr, ctx=self.ctx)
        return self.comms[graph]

    @rpc(NTScalar('?'))
    def create(self, graph):
        return self._get_comm(graph).create()

    @rpc(NTScalar('?'))
    def destroy(self, graph):
        return self._get_comm(graph).destroy()

    @rpc(NTScalar('?'))
    def clear(self, graph):
        return self._get_comm(graph).clear()

    @rpc(NTScalar('?'))
    def reset(self, graph):
        return self._get_comm(graph).reset()

    @rpc(NTScalar('?'))
    def post(self, graph, topic, payload):
        return self._get_comm(graph)._post_dill(topic, dill.loads(payload.tobytes()))

    @rpc(NTScalar('as'))
    def names(self, graph):
        return self._get_comm(graph).names

    @rpc(NTScalar('?'))
    def view(self, graph, name):
        return self._get_comm(graph).view(name)

    @rpc(NTScalar('?'))
    def export(self, graph, name, alias):
        return self._get_comm(graph).export(name, alias)

class PvaExportServer:
    def __init__(self, name, comm_addr, export_addr, aggregate=False):
        self.base = name
        self.ctx = zmq.Context()
        self.export = self.ctx.socket(zmq.SUB)
        self.export.setsockopt_string(zmq.SUBSCRIBE, "")
        self.export.connect(export_addr)
        self.comm = self.ctx.socket(zmq.REQ)
        self.comm.connect(comm_addr)
        self.queue = ThreadedWorkQueue(maxsize=20, workers=1)
        # pva server provider
        self.provider = StaticProvider(name)
        self.rpc_provider = NTURIDispatcher(self.queue,
                                            target=PvaExportRpcHandler(self.ctx, comm_addr),
                                            name="%s:cmd" % self.base,
                                            prefix="%s:cmd:" % self.base)
        self.server_thread = threading.Thread(target=self.server, name='pvaserv')
        self.server_thread.daemon = True
        self.aggregate = aggregate
        self.pvs = {}
        self.ignored = set()
        self.graph_pvbase = "ana"
        self.data_pvbase = "data"
        self.info_pvbase = "info"
        self.cmd_pvs = {'command'}
        self.payload_cmd_pvs = {'add', 'set', 'del'}

    def __enter__(self):
        return self

    def __exit__(self, exc_type, exc_value, traceback):
        self.close()

    def close(self):
        self.ctx.destroy()

    @staticmethod
    def join_pv(*args):
        return ":".join(args)

    def graph_pvname(self, graph, name=None):
        if name is not None:
            return ":".join([self.graph_pvbase, graph, name])
        else:
            return ":".join([self.graph_pvbase, graph])

    def data_pvname(self, graph, name):
        return ":".join([self.graph_pvbase, graph, self.data_pvbase, name])

    def info_pvname(self, name):
        return ":".join([self.info_pvbase, name])

    def find_graph_pvnames(self, graph, names):
        return [name for name in names if name.startswith(self.graph_pvname(graph))]

    def create_pv(self, name, nt, initial, func=None):
        if func is not None:
            pv = SharedPV(nt=nt, initial=initial, handler=PvaExportPutHandler(put=func))
        else:
            pv = SharedPV(nt=nt, initial=initial)
        self.provider.add('%s:%s' % (self.base, name), pv)
        self.pvs[name] = pv

    def create_bytes_pv(self, name, initial, func=None):
        self.create_pv(name, NTBytes(), initial, func=func)

    def valid(self, name, group=None):
        return not name.startswith('_')

    def get_pv_type(self, data):
        if isinstance(data, np.ndarray):
            return NTNDArray()
        elif isinstance(data, bool):
            return NTScalar('?')
        elif isinstance(data, int):
            return NTScalar('l')
        elif isinstance(data, float):
            return NTScalar('d')
        else:
            return NTObject()

    def update_graph(self, graph, data):
        # add the unaggregated version of the pvs
        for key, value in data.items():
            if key in NTGraph.flat_schema:
                name, nttype = NTGraph.flat_schema[key]
                pvname = self.graph_pvname(graph, name)
                if pvname not in self.pvs:
                    self.create_pv(pvname, nttype, value)
                else:
                    self.pvs[pvname].post(value)
        # add the aggregated graph pv if requested
        if self.aggregate:
            pvname = self.graph_pvname(graph)
            if pvname not in self.pvs:
                logger.debug("Creating pv for info on the graph")
                self.create_pv(pvname, NTGraph(), data)
            else:
                self.pvs[pvname].post(data)

    def update_store(self, graph, data):
        # add the unaggregated version of the pvs
        for key, value in data.items():
            if key in NTStore.flat_schema:
                name, nttype = NTStore.flat_schema[key]
                pvname = self.graph_pvname(graph, name)
                if pvname not in self.pvs:
                    self.create_pv(pvname, nttype, value)
                else:
                    self.pvs[pvname].post(value)
        # add the aggregated graph pv if requested
        if self.aggregate:
            pvname = self.graph_pvname(graph, 'store')
            if pvname not in self.pvs:
                logger.debug("Creating pv for info on the store")
                self.create_pv(pvname, NTStore(), data)
            else:
                self.pvs[pvname].post(data)

    def update_heartbeat(self, graph, heartbeat):
        pvname = self.graph_pvname(graph, 'heartbeat')
        if pvname not in self.pvs:
            self.create_pv(pvname, NTScalar('d'), heartbeat.identity)
        else:
            self.pvs[pvname].post(heartbeat.identity)

    def update_info(self, data):
        # add the unaggregated version of the pvs
        for key, value in data.items():
            pvname = self.info_pvname(key)
            if pvname not in self.pvs:
                self.create_pv(pvname, NTScalar('as'), value)
            else:
                self.pvs[pvname].post(value)

    def update_data(self, graph, name, data):
        pvname = self.data_pvname(graph, name)
        if pvname not in self.ignored:
            if pvname not in self.pvs:
                pv_type = self.get_pv_type(data)
                if pv_type is not None:
                    logger.debug("Creating new pv named %s for graph %s", name, graph)
                    self.create_pv(pvname, pv_type, data)
                else:
                    logger.warn("Cannot map type of '%s' from graph '%s' to PV: %s", name, graph, type(data))
                    self.ignored.add(pvname)
            else:
                self.pvs[pvname].post(data)

    def update_destroy(self, graph):
        # close all the pvs associated with the purged graph
        for name in self.find_graph_pvnames(graph, self.pvs):
            logger.debug("Removing pv named %s for graph %s", name, graph)
            self.provider.remove('%s:%s' % (self.base, name))
            del self.pvs[name]
        # remove any ignored pvs associated with the purged graph
        for name in self.find_graph_pvnames(graph, self.ignored):
            self.ignored.remove(name)

    def server(self):
        server = Server(providers=[self.provider, self.rpc_provider])
        with server, self.queue:
            try:
                while True:
                    time.sleep(100)
            except KeyboardInterrupt:
                pass

    def run(self):
        # start the pva server thread
        self.server_thread.start()
        logger.info("Starting PVA data export server")
        while True:
            topic = self.export.recv_string()
            graph = self.export.recv_string()
            exports = self.export.recv_pyobj()
            if topic == 'data':
                for name, data in exports.items():
                    # ignore names starting with '_' - these are private
                    if self.valid(name):
                        self.update_data(graph, name, data)
            elif topic == 'graph':
                self.update_graph(graph, exports)
            elif topic == 'store':
                self.update_store(graph, exports)
            elif topic == 'heartbeat':
                self.update_heartbeat(graph, exports)
            elif topic == 'info':
                self.update_info(exports)
            elif topic == 'destroy':
                self.update_destroy(graph)
            else:
                logger.warn("No handler for topic: %s", topic)


text:
#!/usr/bin/python3

import numpy as np
import pandas as pd

import dask.array as da
import sys
import os

#fluidity_fp  = '/mnt/c/Users/julia/fluidity/fluidity-master/python'
DA_project_fp = '/mnt/c/Users/julia/Documents/Imperial/DA_project'
#sys.path.append(fluidity_fp)
sys.path.append(DA_project_fp)

import vtktools
ug1=vtktools.vtu(DA_project_fp + '/data/LSBU_c_0.vtu')
ug2=vtktools.vtu(DA_project_fp + '/data/LSBU_0.vtu')

print(ug1.GetFieldNames())
print(ug2.GetFieldNames())

#read the values of the tracers and copy in a vector named p
p = ug2.GetVectorField('Velocity')

#create a dask array
p = da.from_array(p, chunks = [5000, 63])
n = len(p)
print(n)

#create background matrix
Background = da.matmul(p, p.T)
Background.compute()
print(Background[0,0])