text:
#!/usr/bin/env python
# coding: utf-8

# # Overview of pMuTT's Core Functionality
# Originally written for Version 1.2.1
#
# Last Updated for Version 1.2.13
#
# ## Topics Covered
#
# - Using constants and converting units using the ``constants`` module
# - Initializing ``StatMech`` objects by specifying all modes and by using ``presets`` dictionary
# - Initializing empirical objects such as ``Nasa`` objects using a ``StatMech`` object or from a previously generated Nasa polynomial
# - Initializing ``Reference`` and ``References`` objects to adjust DFT's reference to more traditional references
# - Input (via Excel) and output ``Nasa`` polynomials to thermdat format
# - Initializing ``Reaction`` objects from strings
#
# ## Useful Links:
#
# - Github: https://github.com/VlachosGroup/pMuTT
# - Documentation: https://vlachosgroup.github.io/pMuTT/index.html
# - Examples: https://vlachosgroup.github.io/pMuTT/examples.html

# ## Constants
# pMuTT has a wide variety of constants to increase readability of the code. See [Constants page][0] in the documentation for supported units.
#
# [0]: https://vlachosgroup.github.io/pmutt/constants.html#constants

# In[1]:

from pmutt import constants as c

1.987
print(c.R('kJ/mol/K'))

print('Some constants')
print('R (J/mol/K) = {}'.format(c.R('J/mol/K')))
print("Avogadro's number = {}\n".format(c.Na))

print('Unit conversions')
print('5 kJ/mol --> {} eV/molecule'.format(c.convert_unit(num=5., initial='kJ/mol', final='eV/molecule')))
print('Frequency of 1000 Hz --> Wavenumber of {} 1/cm\n'.format(c.freq_to_wavenumber(1000.)))

print('See expected inputs, supported units of different constants')
help(c.R)
help(c.convert_unit)

# ## StatMech Objects
# Molecules show translational, vibrational, rotational, electronic, and nuclear modes.
#
#
#
# The [``StatMech``][0] object allows us to specify translational, vibrational, rotational, electronic and nuclear modes independently, which gives flexibility in what behavior you would like.
#
# [0]: https://vlachosgroup.github.io/pmutt/statmech.html#pmutt.statmech.StatMech

# For this example, we will use a butane molecule as an ideal gas:
#
# - translations with no interaction between molecules
# - harmonic vibrations
# - rigid rotor rotations
# - ground state electronic structure -
# - ground state nuclear structure).
#
#

# In[2]:

from ase.build import molecule
from pmutt.statmech import StatMech, trans, vib, rot, elec

butane_atoms = molecule('trans-butane')

'''Translational'''
butane_trans = trans.FreeTrans(n_degrees=3, atoms=butane_atoms)

'''Vibrational'''
butane_vib = vib.HarmonicVib(vib_wavenumbers=[3054.622862, 3047.573455, 3037.53448,
                                              3030.21322, 3029.947329, 2995.758708,
                                              2970.12166, 2968.142985, 2951.122942,
                                              2871.560685, 1491.354921, 1456.480829,
                                              1455.224163, 1429.084081, 1423.153673,
                                              1364.456094, 1349.778994, 1321.137752,
                                              1297.412109, 1276.969173, 1267.783512,
                                              1150.401492, 1027.841298, 1018.203753,
                                              945.310074, 929.15992, 911.661049,
                                              808.685354, 730.986587, 475.287654,
                                              339.164649, 264.682213, 244.584138,
                                              219.956713, 115.923768, 35.56194])

'''Rotational'''
butane_rot = rot.RigidRotor(symmetrynumber=2, atoms=butane_atoms)

'''Electronic'''
butane_elec = elec.GroundStateElec(potentialenergy=-73.7051, spin=0)

'''StatMech Initialization'''
butane_statmech = StatMech(name='butane',
                           trans_model=butane_trans,
                           vib_model=butane_vib,
                           rot_model=butane_rot,
                           elec_model=butane_elec)

H_statmech = butane_statmech.get_H(T=298., units='kJ/mol')
S_statmech = butane_statmech.get_S(T=298., units='J/mol/K')
print('H_butane(T=298) = {:.1f} kJ/mol'.format(H_statmech))
print('S_butane(T=298) = {:.2f} J/mol/K'.format(S_statmech))

# ### Presets
# The [``presets``][0] dictionary stores commonly used models to ease the initialization of [``StatMech``][1] objects. The same water molecule before can be initialized this way instead.
#
# [0]: https://vlachosgroup.github.io/pmutt/statmech.html#presets
# [1]: https://vlachosgroup.github.io/pmutt/statmech.html#pmutt.statmech.StatMech

# In[3]:

from pprint import pprint
from ase.build import molecule
from pmutt.statmech import StatMech, presets

idealgas_defaults = presets['idealgas']
pprint(idealgas_defaults)

# In[4]:

butane_preset = StatMech(name='butane',
                         atoms=molecule('trans-butane'),
                         vib_wavenumbers=[3054.622862, 3047.573455, 3037.53448,
                                          3030.21322, 3029.947329, 2995.758708,
                                          2970.12166, 2968.142985, 2951.122942,
                                          2871.560685, 1491.354921, 1456.480829,
                                          1455.224163, 1429.084081, 1423.153673,
                                          1364.456094, 1349.778994, 1321.137752,
                                          1297.412109, 1276.969173, 1267.783512,
                                          1150.401492, 1027.841298, 1018.203753,
                                          945.310074, 929.15992, 911.661049,
                                          808.685354, 730.986587, 475.287654,
                                          339.164649, 264.682213, 244.584138,
                                          219.956713, 115.923768, 35.56194],
                         symmetrynumber=2,
                         potentialenergy=-73.7051,
                         spin=0,
                         **idealgas_defaults)

H_preset = butane_preset.get_H(T=298., units='kJ/mol')
S_preset = butane_preset.get_S(T=298., units='J/mol/K')
print('H_butane(T=298) = {:.1f} kJ/mol'.format(H_preset))
print('S_butane(T=298) = {:.2f} J/mol/K'.format(S_preset))

# ### Empty Modes
# The [``EmptyMode``][0] is a special object returns 1 for the partition function and 0 for all other thermodynamic properties. This is useful if you do not want any contribution from a mode.
#
# [0]: https://vlachosgroup.github.io/pmutt/statmech.html#empty-mode

# In[5]:

from pmutt.statmech import EmptyMode

empty = EmptyMode()
print('Some EmptyMode properties:')
print('q = {}'.format(empty.get_q()))
print('H/RT = {}'.format(empty.get_HoRT()))
print('S/R = {}'.format(empty.get_SoR()))
print('G/RT = {}'.format(empty.get_GoRT()))

# ## Empirical Objects
# Empirical forms are polynomials that are fit to experimental or *ab-initio data*. These forms are useful because they can be evaluated relatively quickly, so that downstream software is not hindered by evaluation of thermochemical properties.
# However, note that ``StatMech`` objects can calculate more properties than the currently supported empirical objects.

# ### NASA polynomial
# The [``NASA``][0] format is used for our microkinetic modeling software, Chemkin.
#
# #### Initializing Nasa from StatMech
# Below, we initialize the NASA polynomial from the ``StatMech`` object we created earlier.
#
# [0]: https://vlachosgroup.github.io/pmutt/empirical.html#nasa

# In[6]:

from pmutt.empirical.nasa import Nasa

butane_nasa = Nasa.from_model(name='butane',
                              model=butane_statmech,
                              T_low=298.,
                              T_high=800.,
                              elements={'C': 4, 'H': 10},
                              phase='G')

H_nasa = butane_nasa.get_H(T=298., units='kJ/mol')
S_nasa = butane_nasa.get_S(T=298., units='J/mol/K')
print('H_butane(T=298) = {:.1f} kJ/mol'.format(H_nasa))
print('S_butane(T=298) = {:.2f} J/mol/K'.format(S_nasa))

# Although it is not covered here, you can also generate empirical objects from experimental data using the ``.from_data`` method. See [Experimental to Empirical][6] example.
#
# [6]: https://vlachosgroup.github.io/pmutt/examples.html#experimental-to-empirical

# #### Initializing Nasa Directly
# We can also initialize the NASA polynomial if we have the polynomials. Using an entry from the [Reaction Mechanism Generator (RMG) database][0].
#
# [0]: https://rmg.mit.edu/database/thermo/libraries/DFT_QCI_thermo/215/

# In[7]:

import numpy as np

butane_nasa_direct = Nasa(name='butane',
                          T_low=100.,
                          T_mid=1147.61,
                          T_high=5000.,
                          a_low=np.array([ 2.16917742E+00,
                                           3.43905384E-02,
                                          -1.61588593E-06,
                                          -1.30723691E-08,
                                           5.17339469E-12,
                                          -1.72505133E+04,
                                           1.46546944E+01]),
                          a_high=np.array([ 6.78908025E+00,
                                            3.03597365E-02,
                                           -1.21261608E-05,
                                            2.19944009E-09,
                                           -1.50288488E-13,
                                           -1.91058191E+04,
                                           -1.17331911E+01]),
                         elements={'C': 4, 'H': 10},
                         phase='G')

H_nasa_direct = butane_nasa_direct.get_H(T=298., units='kJ/mol')
S_nasa_direct = butane_nasa_direct.get_S(T=298., units='J/mol/K')
print('H_butane(T=298) = {:.1f} kJ/mol'.format(H_nasa_direct))
print('S_butane(T=298) = {:.2f} J/mol/K'.format(S_nasa_direct))

# Compare the results between ``butane_nasa`` and ``butane_nasa_direct`` to the [Wikipedia page for butane][0].
#
# [0]: https://en.wikipedia.org/wiki/Butane_(data_page)

# In[8]:

print('H_nasa = {:.1f} kJ/mol'.format(H_nasa))
print('H_nasa_direct = {:.1f} kJ/mol'.format(H_nasa_direct))
print('H_wiki = -125.6 kJ/mol\n')

print('S_nasa = {:.2f} J/mol/K'.format(S_nasa))
print('S_nasa_direct = {:.2f} J/mol/K'.format(S_nasa_direct))
print('S_wiki = 310.23 J/mol/K')

# Notice the values are very different for ``H_nasa``. This discrepancy is due to:
#
# - different references
# - error in DFT
#
# We can account for this discrepancy by using the [``Reference``][0] and [``References``][1] objects.
#
# [0]: https://vlachosgroup.github.io/pmutt/empirical.html#pmutt.empirical.references.Reference
# [1]: https://vlachosgroup.github.io/pmutt/empirical.html#pmutt.empirical.references.References

# ### Referencing
# To define a reference, you must have:
#
# - enthalpy at some reference temperature (``HoRT_ref`` and ``T_ref``)
# - a ``StatMech`` object
#
# In general, use references that are similar to molecules in your mechanism. Also, the number of reference molecules must be equation to the number of elements (or other descriptor) in the mechanism. [Full description of referencing scheme here][0].
#
# In this example, we will use ethane and propane as references.
#
#
#
# [0]: https://vlachosgroup.github.io/pmutt/referencing.html

# In[9]:

from pmutt.empirical.references import Reference, References

ethane_ref = Reference(name='ethane',
                       elements={'C': 2, 'H': 6},
                       atoms=molecule('C2H6'),
                       vib_wavenumbers=[3050.5296, 3049.8428, 3025.2714,
                                        3024.4304, 2973.5455, 2971.9261,
                                        1455.4203, 1454.9941, 1454.2055,
                                        1453.7038, 1372.4786, 1358.3593,
                                        1176.4512, 1175.507, 992.55,
                                        803.082, 801.4536, 298.4712],
                       symmetrynumber=6,
                       potentialenergy=-40.5194,
                       spin=0,
                       T_ref=298.15,
                       HoRT_ref=-33.7596,
                       **idealgas_defaults)

propane_ref = Reference(name='propane',
                        elements={'C': 3, 'H': 8},
                        atoms=molecule('C3H8'),
                        vib_wavenumbers=[3040.9733, 3038.992, 3036.8071,
                                         3027.6062, 2984.8436, 2966.1692,
                                         2963.4684, 2959.7431, 1462.5683,
                                         1457.4221, 1446.858, 1442.0357,
                                         1438.7871, 1369.6901, 1352.6287,
                                         1316.215, 1273.9426, 1170.4456,
                                         1140.9699, 1049.3866, 902.8507,
                                         885.3209, 865.5958, 735.1924,
                                         362.7372, 266.3928, 221.4547],
                        symmetrynumber=2,
                        potentialenergy=-57.0864,
                        spin=0,
                        T_ref=298.15,
                        HoRT_ref=-42.2333,
                       **idealgas_defaults)

refs = References(references=[ethane_ref, propane_ref])
print(refs.offset)

# Passing the ``References`` object when we make our ``Nasa`` object produces a value closer to the one listed above.

# In[10]:

butane_nasa_ref = Nasa.from_model(name='butane',
                                  model=butane_statmech,
                                  T_low=298.,
                                  T_high=800.,
                                  elements={'C': 4, 'H': 10},
                                  references=refs)
H_nasa_ref = butane_nasa_ref.get_H(T=298., units='kJ/mol')
S_nasa_ref = butane_nasa_ref.get_S(T=298., units='J/mol/K')
print('H_butane(T=298) = {:.1f} kJ/mol'.format(H_nasa_ref))
print('S_butane(T=298) = {:.2f} J/mol/K'.format(S_nasa_ref))

# ## Input and Output
# ### Excel
# Encoding each object in Python can be tedious and so you can read from Excel spreadsheets using [``pmutt.io.excel.read_excel``][0]. Note that this function returns a list of dictionaries. This output allows you to initialize whichever object you want. There are also special rules that depend on the header name.
#
# [0]: https://vlachosgroup.github.io/pmutt/io.html?highlight=read_excel#pmutt.io.excel.read_excel

# In[11]:

import os
from pathlib import Path
from pmutt.io.excel import read_excel

# Find the location of Jupyter notebook
# Note that normally Python scripts have a __file__ variable but Jupyter notebook doesn't.
# Using pathlib can overcome this limiation
try:
    notebook_folder = os.path.dirname(__file__)
except NameError:
    notebook_folder = Path().resolve()
os.chdir(notebook_folder)

# The Excel spreadsheet is located in the same folder as the Jupyter notebook
refs_path = os.path.join(notebook_folder, 'refs.xlsx')

refs_data = read_excel(refs_path)
pprint(refs_data)

# Initialize using \*\*kwargs syntax.

# In[12]:

ref_list = []
for record in refs_data:
    ref_list.append(Reference(**record))

refs_excel = References(references=ref_list)
print(refs_excel.offset)

# Butane can be initialized in a similar way.

# In[13]:

# The Excel spreadsheet is located in the same folder as the Jupyter notebook
butane_path = os.path.join(notebook_folder, 'butane.xlsx')
butane_data = read_excel(butane_path)[0] # [0] accesses the butane data

butane_excel = Nasa.from_model(T_low=298.,
                               T_high=800.,
                               references=refs_excel,
                               **butane_data)

H_excel = butane_excel.get_H(T=298., units='kJ/mol')
S_excel = butane_excel.get_S(T=298., units='J/mol/K')
print('H_butane(T=298) = {:.1f} kJ/mol'.format(H_excel))
print('S_butane(T=298) = {:.2f} J/mol/K'.format(S_excel))

# ### Thermdat
# The thermdat format uses NASA polynomials to represent several species. It has a very particular format so doing it manually is error-prone. You can write a list of ``Nasa`` objects to thermdat format using [``pmutt.io.thermdat.write_thermdat``][0].
#
# [0]: https://vlachosgroup.github.io/pmutt/io.html#pmutt.io.thermdat.write_thermdat

# In[14]:

from pmutt.io.thermdat import write_thermdat

# Make Nasa objects from previously defined ethane and propane
ethane_nasa = Nasa.from_model(name='ethane',
                              phase='G',
                              T_low=298.,
                              T_high=800.,
                              model=ethane_ref.model,
                              elements=ethane_ref.elements,
                              references=refs)
propane_nasa = Nasa.from_model(name='propane',
                               phase='G',
                               T_low=298.,
                               T_high=800.,
                               model=propane_ref.model,
                               elements=propane_ref.elements,
                               references=refs)
nasa_species = [ethane_nasa, propane_nasa, butane_nasa]

# Determine the output path and write the thermdat file
thermdat_path = os.path.join(notebook_folder, 'thermdat')
write_thermdat(filename=thermdat_path, nasa_species=nasa_species)

THERMO ALL
       100       500      1500
ethane          20190122C   2H   6          G298.0     800.0     502.9         1
-6.84729317E-01 2.54150584E-02-1.02900193E-05-1.15304090E-09 1.73800327E-12    2
-1.08866341E+04 2.42655635E+01 6.77655563E+00-3.31037518E-02 1.63456769E-04    3
-2.32569163E-07 1.18361103E-10-1.16548936E+04-6.74224205E+00                   4
propane         20190122C   3H   8          G298.0     800.0     502.9         1
-2.54716915E+00 4.36640749E-02-2.65978975E-05 6.89956646E-09 7.68567634E-14    2
-1.35353831E+04 3.50227651E+01 8.09253651E+00-4.01657054E-02 2.23439464E-04    3
-3.27632587E-07 1.69405578E-10-1.46259783E+04-9.14554121E+00                   4
butane          20190122C   4H  10          G298.0     800.0     502.9         1
-1.98482876E+00 5.99910769E-02-4.67983727E-05 2.14853318E-08-4.31005098E-12    2
-8.15415828E+05 3.48691129E+01 9.98107805E+00-3.50556935E-02 2.38938732E-04    3
-3.63709779E-07 1.92062158E-10-8.16632302E+05-1.47065846E+01                   4
END
# Similarly, [``pmutt.io.thermdat.read_thermdat``][0] reads thermdat files.
#
# [0]: https://vlachosgroup.github.io/pmutt/io.html#pmutt.io.thermdat.read_thermdat

# ## Reactions
# You can also evaluate reactions properties. The most straightforward way to do this is to initialize using strings.

# In[15]:

from pmutt.io.thermdat import read_thermdat
from pmutt import pmutt_list_to_dict
from pmutt.reaction import Reaction

# Get a dictionary of species
thermdat_H2O_path = os.path.join(notebook_folder, 'thermdat_H2O')
species_list = read_thermdat(thermdat_H2O_path)
species_dict = pmutt_list_to_dict(species_list)

# Initialize the reaction
rxn_H2O = Reaction.from_string('H2 + 0.5O2 = H2O', species=species_dict)

# Calculate reaction properties
H_rxn = rxn_H2O.get_delta_H(T=298., units='kJ/mol')
S_rxn = rxn_H2O.get_delta_S(T=298., units='J/mol/K')
print('H_rxn(T=298) = {:.1f} kJ/mol'.format(H_rxn))
print('S_rxn(T=298) = {:.2f} J/mol/K'.format(S_rxn))

# ## Exercise
# Write a script to calculate the Enthalpy of adsorption (in kcal/mol) of H2O on Cu(111) at T = 298 K. Some important details are given below.
#
# ### Information Required
# #### H2O:
# - ideal gas
# - atoms: You can use "ase.build.molecule" to generate a water molecule like we did with ethane, propane, and butane.
# - vibrational wavenumbers (1/cm): 3825.434, 3710.2642, 1582.432
# - potential energy (eV): -14.22393533
# - spin: 0
# - symmetry number: 2
#
# #### Cu(111):
# - only electronic modes
# - potential energy (eV): -224.13045381
#
# #### H2O+Cu(111):
# - electronic and harmonic vibration modes
# - potential energy (eV): -238.4713854
# - vibrational wavenumbers (1/cm): 3797.255519, 3658.895695, 1530.600295, 266.366130, 138.907356, 63.899768, 59.150454, 51.256019, -327.384554 (negative numbers represent imaginary frequencies. The default behavior of pMuTT is to ignore these frequencies when calculating any thermodynamic property)
#
# #### Reaction:
# H2O + Cu(111) --> H2O+Cu(111)
#
# ### Solution:

# In[16]:

from ase.build import molecule
from pmutt.statmech import StatMech, presets
from pmutt.reaction import Reaction

# Using dictionary since later I will initialize the reaction with a string
species = {
    'H2O(g)': StatMech(atoms=molecule('H2O'),
                       vib_wavenumbers=[3825.434, 3710.2642, 1582.432],
                       potentialenergy=-14.22393533,
                       spin=0,
                       symmetrynumber=2,
                       **presets['idealgas']),
    '*': StatMech(potentialenergy=-224.13045381,
                  **presets['electronic']),
    'H2O*': StatMech(potentialenergy=-238.4713854,
                     vib_wavenumbers=[3797.255519,
                                      3658.895695,
                                      1530.600295,
                                      266.366130,
                                      138.907356,
                                      63.899768,
                                      59.150454,
                                      51.256019,
                                      -327.384554], #Imaginary frequency!
                     **presets['harmonic']),
}

rxn = Reaction.from_string('H2O(g) + * = H2O*', species)
del_H = rxn.get_delta_H(T=298., units='kcal/mol')
print('del_H = {:.2f} kcal/mol'.format(del_H))

# In[ ]:

max_stars_repo_name: wittregr/pMuTT

max_issues_repo_name: wittregr/pMuTT

max_forks_repo_name: wittregr/pMuTT


text:
from PIL import Image
from math import floor
import numpy as np
from sys import argv

def main():

	im = Image.open('imagen_Karina.jpg')
    	ancho,alto = im.size
    	imagen = im.load()
    	im.show()

	print "Tamano original: ", ancho,alto
	x = int(argv[1])
	y = int(argv[2])

    	mide_ancho = (1.0*x)/ancho
    	mide_alto = (1.0*y)/alto

    	n_ancho = int(ancho*mide_ancho)
    	n_alto = int(alto*mide_alto)

    	nueva = Image.new('RGB',(n_ancho,n_alto))

	print "Nuevo tamano: ", n_ancho,n_alto
    	otra = nueva.load()

    	res = aum_dis(nueva,imagen,otra,mide_ancho,mide_alto,n_ancho,n_alto)
    	res.show()

def aum_dis(im,imagen,otra,mide_ancho,mide_alto,n_ancho,n_alto):

    for a in range(n_ancho):
        for b in range(n_alto):
            otra[a,b] = imagen[int(a/mide_ancho),int(b/mide_alto)]

    return im

if __name__ == '__main__':
    main()

ext: py

reason: import numpy


text:
import base
import numpy as np
import tensorflow as tf
import copy

### Basic modules describing the core search and computational functionality.
class BasicModule(object):
    """Basic search and compilation functionality for basic modules.

    Basic modules are composed of parameters and hyperparameters. They are not
    composite in the sense that do not have submodules as it is the case of
    Concat and Or modules. This class is not meant to be used directly.
    It has to be inherited and the functions that compute the output dimension
    and compile the module to tensorflow has to be implemented for the specific
    module considered. See Affine and Dropout class definitions for examples
    of how to extend this class.

    """

    def __init__(self):
        self.order = []
        self.domains = []
        self.chosen = []
        self.in_d = None
        self.scope = None
        self.namespace_id = None

    # for search
    def initialize(self, in_d, scope):
        if self.in_d != None or self.scope != None:
            raise AssertionError
        #print(in_d)
        if len(in_d) < 1:
            raise ValueError

        self.in_d = in_d
        self.scope = scope

        # registers itself in the scope
        prefix = self.__class__.__name__
        name = scope.get_valid_name(prefix)
        self.namespace_id = name
        scope.register_namespace(name)
        propagate(self)

    def get_outdim(self):
        raise NotImplemented

    def is_specified(self):
        return (self.scope != None and
                self.in_d != None and
                len(self.chosen) == len(self.order))

    def get_choices(self):
        assert not self.is_specified()

        i = len(self.chosen)
        return (self.order[i], self.domains[i])

    def choose(self, choice_i):
        assert not self.is_specified()

        # print self.chosen, choice_i
        if len(self.domains[len(self.chosen)]) <= choice_i:
            raise ValueError

        self.chosen.append(choice_i)
        propagate(self)

    # for printing
    def repr_program(self):
        name = self.__class__.__name__
        return (name, ) + tuple(self.domains)

    def repr_model(self):
        name = self.__class__.__name__
        vals = [dm[i] for i, dm in zip(self.chosen, self.domains)]
        r = (name, ) + tuple(vals)
        return (r, )

    # for generating the tensorflow code
    def compile(self, in_x, train_feed, eval_feed):
        raise NotImplemented

class Empty(BasicModule):
    """Empty Module.

    Compiles to a wire. Directly passes the input to the output without
    any transformation.

    """

    def get_outdim(self):
        return self.in_d

    def compile(self, in_x, train_feed, eval_feed):
        return in_x

class Dropout(BasicModule):
    """Dropout module.

    Dropout has different training behavior depending on whether the network
    is being run on a training or on evaluation phase. Getting the desired
    behavior in each case is achieved through the use of scopes.

    """

    def __init__(self, ps):
        super(Dropout, self).__init__()
        self.order.append("keep_prob")
        self.domains.append(ps)

    def get_outdim(self):
        return self.in_d

    def compile(self, in_x, train_feed, eval_feed):
        p_name = self.namespace_id + '_' + self.order[0]
        p_var = tf.placeholder(tf.float32, name=p_name)

        # during training the value of the dropout probability (keep_prob) is
        # set to the actual chosen value.
        # during evalution, it is set to 1.0.
        p_val = self.domains[0][self.chosen[0]]
        train_feed[p_var] = p_val
        eval_feed[p_var] = 1.0
        out_y = tf.nn.dropout(in_x, p_var)
        return out_y

class BatchNormalization(BasicModule):
    def get_outdim(self):
        return self.in_d

    def compile(self, in_x, train_feed, eval_feed):
        p_name = self.namespace_id + '_' + "IsTraining"
        p_var = tf.placeholder(tf.bool, name=p_name)
        train_feed[p_var] = True
        eval_feed[p_var] = False
        out_y = tf.contrib.layers.batch_norm(
                    in_x, decay=0.9, is_training=p_var, updates_collections=None)
        return out_y

class ReLU(BasicModule):
    def get_outdim(self):
        return self.in_d

    def compile(self, in_x, train_feed, eval_feed):
        out_y = tf.nn.relu(in_x)
        return out_y

class Sigmoid(BasicModule):
    def get_outdim(self):
        return self.in_d

    def compile(self, in_x, train_feed, eval_feed):
        out_y = tf.nn.sigmoid(in_x)
        return out_y

class Tanh(BasicModule):
    def get_outdim(self):
        return self.in_d

    def compile(self, in_x, train_feed, eval_feed):
        out_y = tf.nn.tanh(in_x)
        return out_y

class Affine(BasicModule):
    def __init__(self, ms, param_init_fns):
        super(Affine, self).__init__()
        self.order.extend(["num_hidden_units", "param_init_fn"])
        self.domains.extend([ms, param_init_fns])

    def get_outdim(self):
        return (self.domains[0][self.chosen[0]], )

    def compile(self, in_x, train_feed, eval_feed):
        n = np.product(self.in_d)
        m, param_init_fn = [dom[i] for (dom, i) in zip(self.domains, self.chosen)]

        #sc = np.sqrt(6.0) / np.sqrt(m + n)
        #W = tf.Variable(tf.random_uniform([n, m], -sc, sc))
        W = tf.Variable( param_init_fn( [n, m] ) )
        b = tf.Variable(tf.zeros([m]))

        # if the number of input dimensions is larger than one, flatten the
        # input and apply the affine transformation.
        if len(self.in_d) > 1:
            in_x_flat = tf.reshape(in_x, shape=[-1, n])
            out_y = tf.add(tf.matmul(in_x_flat, W), b)
        else:
            out_y = tf.add(tf.matmul(in_x, W), b)
        return out_y

# computes the output dimension based on the padding scheme used.
# this comes from the tensorflow documentation
def compute_padded_dims(in_height, in_width, window_len, stride, padding):

    if padding == "SAME":
        out_height = int(np.ceil(float(in_height) / float(stride)))
        out_width  = int(np.ceil(float(in_width) / float(stride)))
    elif padding == "VALID":
        out_height = int(np.ceil(float(in_height - window_len + 1) / float(stride)))
        out_width  = int(np.ceil(float(in_width - window_len + 1) / float(stride)))
    else:
        raise ValueError

    return (out_height, out_width)

class MaxPooling2D(BasicModule):
    def __init__(self, window_lens, strides, paddings):
        super(MaxPooling2D, self).__init__()

        self.order.extend(["window_len", "stride", "padding"])
        self.domains.extend([window_lens, strides, paddings])

    # does additional error checking on the dimension.
    def initialize(self, in_d, scope):
        if len(in_d) != 3:
            raise ValueError
        else:
            super(MaxPooling2D, self).initialize(in_d, scope)

    def get_outdim(self):
        in_height, in_width, in_nchannels = self.in_d
        window_len, stride, padding = [dom[i]
                for (dom, i) in zip(self.domains, self.chosen)]

        out_height, out_width = compute_padded_dims(
            in_height, in_width, window_len, stride, padding)
        out_d = (out_height, out_width, in_nchannels)

        return out_d

    def compile(self, in_x, train_feed, eval_feed):
        in_height, in_width, in_nchannels = self.in_d
        window_len, stride, padding = [dom[i]
                for (dom, i) in zip(self.domains, self.chosen)]

        out_y = tf.nn.max_pool(
            in_x, ksize=[1, window_len, window_len, 1],
            strides=[1, stride, stride, 1], padding=padding)

        return out_y

# NOTE: for now, this is close to a replication of the max pooling layer.
# this may change later if we can capture most of the pooling layers in the
# same format. for now, an auxiliary function is provided.
class AvgPooling2D(BasicModule):
    def __init__(self, window_lens, strides, paddings):
        super(AvgPooling2D, self).__init__()

        self.order.extend(["window_len", "stride", "padding"])
        self.domains.extend([window_lens, strides, paddings])

    # does additional error checking on the dimension.
    def initialize(self, in_d, scope):
        if len(in_d) != 3:
            raise ValueError
        else:
            super(AvgPooling2D, self).initialize(in_d, scope)

    def get_outdim(self):
        in_height, in_width, in_nchannels = self.in_d
        window_len, stride, padding = [dom[i]
                for (dom, i) in zip(self.domains, self.chosen)]

        out_height, out_width = compute_padded_dims(
            in_height, in_width, window_len, stride, padding)
        out_d = (out_height, out_width, in_nchannels)

        return out_d

    def compile(self, in_x, train_feed, eval_feed):
        in_height, in_width, in_nchannels = self.in_d
        window_len, stride, padding = [dom[i]
                for (dom, i) in zip(self.domains, self.chosen)]

        out_y = tf.nn.avg_pool(
            in_x, ksize=[1, window_len, window_len, 1],
            strides=[1, stride, stride, 1], padding=padding)

        return out_y

class Conv2D(BasicModule):

    def __init__(self, filter_numbers, filter_lens, strides, paddings,
            param_init_fns):
        super(Conv2D, self).__init__()

        self.order.extend(["filter_number", "filter_len", "stride", "padding",
            "param_init_fn"])
        self.domains.extend([filter_numbers, filter_lens, strides, paddings,
            param_init_fns])

    # does additional error checking on the dimension.
    def initialize(self, in_d, scope):
        #print(in_d)
        if len(in_d) != 3:
            raise ValueError
        else:
            super(Conv2D, self).initialize(in_d, scope)

    def get_outdim(self):
        in_height, in_width, in_nchannels = self.in_d
        nfilters, filter_len, stride, padding, _ = [dom[i]
                for (dom, i) in zip(self.domains, self.chosen)]

        out_height, out_width = compute_padded_dims(
                in_height, in_width, filter_len, stride, padding)
        out_d = (out_height, out_width, nfilters)

        return out_d

    def compile(self, in_x, train_feed, eval_feed):
        in_height, in_width, in_nchannels = self.in_d
        nfilters, filter_len, stride, padding, param_init_fn = [dom[i]
                for (dom, i) in zip(self.domains, self.chosen)]

        # Creation and initialization of the parameters. Should take size of
        # the filter into account.
        W = tf.Variable(
                param_init_fn( [filter_len, filter_len, in_nchannels, nfilters]) )
        b = tf.Variable(tf.zeros([nfilters]))

        # create the output and add the bias.
        out_yaux = tf.nn.conv2d(in_x, W, strides=[1, stride, stride, 1], padding=padding)
        out_y = tf.nn.bias_add(out_yaux, b)

        #print(in_x.get_shape(), self.get_outdim(), out_y.get_shape())

        return out_y

class UserHyperparams(BasicModule):
    """Used by the user to specify an additional set of hyperparameters that
    the user also wants to search over. The behavior of the program with
    respect to these hyperparameters is determined by the user.
    """
    def __init__(self, order, domains):
        super(UserHyperparams, self).__init__()
        self.order = order
        self.domains = domains

    def get_outdim(self):
        return self.in_d

    # only publishes the values of the hyperparameters to the scope.
    def compile(self, in_x, train_feed, eval_feed):
        namespace = self.scope.get_namespace(self.namespace_id)
        namespace["hyperp_names"] = self.order
        namespace["choices"] = tuple(self.chosen)
        namespace["hyperp_vals"] = [dom[i] for (dom, i) in zip(self.domains, self.chosen)]

        return in_x

### Auxiliary function for dimension and choice propagation.
def propagate_seq(bs, i):
    """ Propagates choices in a sequence of modules.

    If the module in the current position of the sequence is specified, we can
    initialize the next module in the sequence (if there is any), and go to
    the next module in the chain if the initialized module becomes specified.

    """

    while bs[i].is_specified():
        prev_scope = bs[i].scope
        prev_out_d = bs[i].get_outdim()

        i += 1
        if i < len(bs):
            bs[i].initialize(prev_out_d, prev_scope)
        else:
            break
    return i

def propagate(b):
    """ Propagates choices in a module.

    While the module is in a state where there is only one option available for
    the next choice, we take that choice. This function leaves the module
    specified or in a state where there multiple choices. This function will
    typically be called by the submodule when initialize or choose is called
    on that module.

    """
    while not b.is_specified() and len(b.get_choices()[1]) == 1:
        b.choose(0)

### Composite modules that take other modules as input.
class Concat:
    def __init__(self, bs):
        if len(bs) == 0:
            raise ValueError

        self.bs = bs
        self.in_d = None
        self.scope = None

    def initialize(self, in_d, scope):
        self.in_d = in_d
        self.scope = scope

        self.b_index = 0
        self.bs[0].initialize(self.in_d, self.scope)
        self.b_index = propagate_seq(self.bs, self.b_index)

    def get_outdim(self):
        return self.bs[-1].get_outdim()

    def is_specified(self):
        return self.bs[-1].is_specified()

    def get_choices(self):
        return self.bs[self.b_index].get_choices()

    def choose(self, choice_i):
        self.bs[self.b_index].choose(choice_i)
        self.b_index = propagate_seq(self.bs, self.b_index)

    def repr_program(self):
        name = self.__class__.__name__
        args = [b.repr_program() for b in self.bs]
        return (name, ) + tuple(args)

    def repr_model(self):
        vals = []
        for b in self.bs:
            vals.extend(b.repr_model())
        return tuple(vals)

    def compile(self, in_x, train_feed, eval_feed):
        out_y_prev = in_x
        for b in self.bs:
            out_y_prev = b.compile(out_y_prev, train_feed, eval_feed)
        return out_y_prev

class Or:
    def __init__(self, bs):
        if len(bs) == 0:
            raise ValueError

        self.order = ["or_branch"]
        self.domains = [ range(len(bs)) ]
        self.bs = bs
        self.in_d = None
        self.chosen = []
        self.scope = None

    def initialize(self, in_d, scope):
        self.in_d = in_d
        self.scope = scope
        propagate(self)

    def get_outdim(self):
        i = self.chosen[0]

        return self.bs[i].get_outdim()

    def is_specified(self):
        if len(self.chosen) == 1:
            i = self.chosen[0]
            return self.bs[i].is_specified()
        else:
            return False

    def get_choices(self):
        if len(self.chosen) == 0:
            return (self.order[0], self.domains[0])
        else:
            i = self.chosen[0]
            return self.bs[i].get_choices()

    def choose(self, choice_i):
        if len(self.chosen) == 0:
            self.chosen.append(choice_i)
            self.bs[choice_i].initialize(self.in_d, self.scope)
        else:
            b_i = self.chosen[0]
            self.bs[b_i].choose(choice_i)

    def repr_program(self):
        name = self.__class__.__name__
        args = [b.repr_program() for b in self.bs]
        return (name, ) + tuple(args)

    def repr_model(self):
        vals = self.bs[self.chosen[0]].repr_model()
        return vals

    def compile(self, in_x, train_feed, eval_feed):
        out_y = self.bs[self.chosen[0]].compile(in_x, train_feed, eval_feed)
        return out_y

class Repeat:
    """Repeat module.

    Takes a module as input and repeats it some number of times. The number of
    repeats is itself an hyperparameter. Hyperparameters of each repeat are not
    tied across repeats. See the RepeatTied module for tied hyperparameters
    across repeats.

    """
    def __init__(self, b, ks):
        if any([k < 1 for k in ks]):
            raise ValueError

        self.order = ["num_repeats"]
        self.domains = [ks]
        self.b = b

        # used during search
        self.in_d = None
        self.chosen = []
        self.active_bs = None
        self.b_index = None

    def initialize(self, in_d, scope):
        self.in_d = in_d
        self.scope = scope
        propagate(self)

    def get_outdim(self):
        return self.active_bs[-1].get_outdim()

    def is_specified(self):
        if len(self.chosen) == 1:
            return self.active_bs[-1].is_specified()
        else:
            return False

    def get_choices(self):
        if len(self.chosen) == 0:
            return (self.order[0], self.domains[0])
        else:
            return self.active_bs[self.b_index].get_choices()

    def choose(self, choice_i):
        if len(self.chosen) == 0:
            self.chosen.append(choice_i)
            k = self.domains[0][choice_i]

            self.b_index = 0
            self.active_bs = [copy.deepcopy(self.b) for _ in xrange(k)]
            self.active_bs[0].initialize(self.in_d, self.scope)
        else:
            self.active_bs[self.b_index].choose(choice_i)

        self.b_index = propagate_seq(self.active_bs, self.b_index)

    def repr_program(self):
        name = self.__class__.__name__
        args = self.b.repr_program()

        return (name, ) + tuple(args) + tuple(self.domains)

    def repr_model(self):
        vals = []
        for b in self.active_bs:
            vals.extend(b.repr_model())
        return tuple(vals)

    def compile(self, in_x, train_feed, eval_feed):
        out_y_prev = in_x
        for b in self.active_bs:
            out_y_prev = b.compile(out_y_prev, train_feed, eval_feed)
        return out_y_prev

class RepeatTied:
    """RepeatTied module.

    Like a Repeat module, but now the hyperparameters are tied across repeats.
    Note that only the hyperparameters, not the parameters, are tied.
    Tying parameters is done through other route.

    """
    def __init__(self, b, ks):
        if any([k < 1 for k in ks]):
            raise ValueError

        self.order = ["num_repeats"]
        self.domains = [ks]
        self.b = b

        # used during search
        self.in_d = None
        self.chosen = []
        self.active_bs = None
        self.b_index = None
        # this field is for the reuse in setting the repeats b
        self.b0_choose_hist = []

    def initialize(self, in_d, scope):
        self.in_d = in_d
        self.scope = scope
        propagate(self)

    def get_outdim(self):
        return self.active_bs[-1].get_outdim()

    def is_specified(self):
        if len(self.chosen) == 1:
            return self.active_bs[-1].is_specified()
        else:
            return False

    def get_choices(self):
        if len(self.chosen) == 0:
            return (self.order[0], self.domains[0])
        else:
            return self.active_bs[self.b_index].get_choices()

    def choose(self, choice_i):
        if len(self.chosen) == 0:
            self.chosen.append(choice_i)
            k = self.domains[0][choice_i]

            self.b_index = 0
            self.active_bs = [copy.deepcopy(self.b) for _ in xrange(k)]
            self.active_bs[0].initialize(self.in_d, self.scope)
        else:
            self.active_bs[self.b_index].choose(choice_i)
            self.b0_choose_hist.append(choice_i)

        self.b_index = propagate_seq(self.active_bs, self.b_index)

        # as soon as the pointer moves to 1, set all others by copying the
        # choices.
        if self.b_index == 1:
            for b in self.active_bs[1:]:
                for ch in self.b0_choose_hist:
                    b.choose(ch)

                # this should advance the pointer by one.
                self.b_index = propagate_seq(self.active_bs, self.b_index)

            # it should be specified at the end of this.
            assert self.is_specified()

    def repr_program(self):
        name = self.__class__.__name__
        args = self.b.repr_program()

        return (name, ) + tuple(args) + tuple(self.domains)

    def repr_model(self):
        vals = []
        for b in self.active_bs:
            vals.extend(b.repr_model())
        return tuple(vals)

    def compile(self, in_x, train_feed, eval_feed):
        out_y_prev = in_x
        for b in self.active_bs:
            out_y_prev = b.compile(out_y_prev, train_feed, eval_feed)
        return out_y_prev

class Residual:
    """ Residual skip connection.

    Introduces a skip connection between the module passed as argument.
    The module taken as argument can have hyperparameters to be specified.

    If the input and output do not have the same dimensions, padding needs
    to be done for the results to be combined in a sum or product. We briefly
    discuss the different cases:
    (1)| both input and output have the same number of dimensions and the same
    sizes for paired dimensions.
    => simply do the entrywise operation without doing any changes.
    (2)| both input and output have the same number of dimensions, but they
    have different sizes for paired dimensions.
    => pad the smallest dimensions on either input or output such that the
    result after padding can be combined; both input and output can be changed
    in this case if none strictly dominates the other in terms of dimensions.
    (3)| input and ouput have different number of dimensions.
    => perhaps the most straightforward solution is to flatten both input and
    output and combine the flattened versions. another possibility is to add
    extra dimensions and pad the smallest dimensions with zeros.

    The most straightforward solutions have been implemented for now.

    """

    def __init__(self, b):
        self.b = b
        self.in_d = None
        self.scope = None

    def initialize(self, in_d, scope):
        self.in_d = in_d
        self.scope = scope
        self.b.initialize(in_d, scope)
        propagate(self.b)

    def get_outdim(self):
        #assert in_x == self.b.get_outdim()
        # relaxing input dimension equal to output dimension. taking into
        # account the padding scheme considered.
        out_d_b = self.b.get_outdim()
        in_d = self.in_d

        if len(out_d_b) == len(in_d):
            out_d = tuple(
                [max(od_i, id_i) for (od_i, id_i) in zip(out_d_b, in_d)])

        else:
            # flattens both input and output.
            out_d_b_flat = np.product(out_d_b)
            in_d_flat = np.product(in_d)
            out_d = (max(out_d_b_flat, in_d_flat) ,)

        return out_d

    def is_specified(self):
        return self.b.is_specified()

    def get_choices(self):
        return self.b.get_choices()

    def choose(self, choice_i):
        self.b.choose(choice_i)

    def repr_program(self):
        name = self.__class__.__name__
        args = self.b.repr_program()

        return (name, args)

    def repr_model(self):
        name = self.__class__.__name__
        b = self.b.repr_model()
        r = (name, b[0])

        return (r, )

    def compile(self, in_x, train_feed, eval_feed):
        # NOTE: this function requires that target dims dominate (bigger or
        # equal component wise) in_dims. this is the case in how it is used
        # currently in the code below.
        compute_padding_fn = lambda in_dims, out_dims: [ [0, max(0, od_i - id_i)]
                for (id_i, od_i) in zip(in_dims, out_dims) ]

        out_d_b = self.b.get_outdim()
        in_d = self.in_d
        out_d = self.get_outdim()
        out_y_b = self.b.compile(in_x, train_feed, eval_feed)

        # if the number of dimensions is not the same, flatten both.
        if len(out_d) != len(in_d):
            out_d_b_flat = np.product(out_d_b)
            out_y_b = tf.reshape(out_y_b, [-1, out_d_b_flat])
            out_d_b = (out_d_b_flat, )

            in_d_flat = np.product(in_d)
            in_x = tf.reshape(in_x, [-1,  in_d_flat])
            in_d = (in_d_flat, )

        # computing the padding for both b and the input.
        # NOTE: adds no padding to the data dimension (i.e., the initial [0, 0])
        paddings_b = [[0, 0]] + compute_padding_fn(out_d_b, out_d)
        out_y_b_padded = tf.pad(out_y_b, paddings_b, "CONSTANT")

        paddings_in = [[0, 0]] + compute_padding_fn(in_d, out_d)
        in_x_padded = tf.pad(in_x, paddings_in, "CONSTANT")

        # finally combine the results with the appropriate dimensions.
        out_y = out_y_b_padded + in_x_padded

        return out_y

class Squeeze(BasicModule):
    def __init__(self):
        super(Squeeze, self).__init__()

    def initialize(self, in_d, scope):
        super(Squeeze, self).initialize(in_d, scope)

    def get_outdim(self):
        out_d = filter(lambda a: a != 1, self.in_d)
        return out_d

    def compile(self, in_x, train_feed, eval_feed):
        out_y = tf.squeeze(in_x)
        return out_y

# it is kind of like the previous one.
class ChoiceBisection:
    """Does bissection on all the hyperparameters of module taken as argument.
    This can be useful to increase sharing between hyperparameter values in
    approaches such as MCTS.
    """
    def __init__(self, b):
        self.b = b
        self.is_bisecting = False
        self.left = None
        self.right = None
        self.hist = []
        self.cur_name = None
        self.cur_vals = None
        # hist is kept to create the create a name for the binary choice with
        # a suffix that corresponds to the sequence of binary decisions done
        # so far.

    def initialize(self, in_d, scope):
        self.b.initialize(in_d, scope)

    def get_outdim(self):
        return self.b.get_outdim()

    def is_specified(self):
        return self.b.is_specified()

    def get_choices(self):
        # only potentially called in the first call to get_choices.
        if not self.is_bisecting:
            self.cur_name, self.cur_vals = self.b.get_choices()
            self.left = 0
            self.right = len(self.cur_vals)
            self.is_bisecting = True
            assert len(self.cur_vals) > 1

        name = self.cur_name + '_' + "".join(map(str, self.hist))
        choices = (0, 1)

        return (name, choices)

    def choose(self, choice_i):
        if not self.is_bisecting:
            self.cur_name, self.cur_vals = self.b.get_choices()
            self.left = 0
            self.right = len(self.cur_vals)
            self.is_bisecting = True
            assert len(self.cur_vals) > 1

        assert self.is_bisecting and choice_i == 0 or choice_i == 1

        # doing the bissection
        mid = int( (self.left + self.right) / 2.0 )
        if choice_i == 0:
            self.right = mid
        else:
            self.left = mid

        self.hist.append(choice_i)

        # comitting to the choice when there is a single one left.
        if self.right - self.left == 1:
            self.b.choose(self.left)
            self.is_bisecting = False
            self.left = None
            self.right = None
            self.hist = []

    def repr_program(self):
        name = self.__class__.__name__
        args = self.b.repr_program()

        return (name, args)

    def repr_model(self):
        return self.b.repr_model()

    # for generating the tensorflow code
    def compile(self, in_x, train_feed, eval_feed):
        return self.b.compile(in_x, train_feed, eval_feed)

# auxiliary functions for some common operations.
def Optional_fn(b):
    return Or([Empty(), b])

def Nonlinearity_fn(nonlin_types):
    bs = []

    for t in nonlin_types:
        if t == "relu":
            b = ReLU()
        elif t == 'sigmoid':
            b = Sigmoid()
        elif t == 'tanh':
            b = Tanh()
        else:
            raise ValueError

        bs.append(b)

    return Or(bs)

def Pooling2D_fn(pooling_types, window_lens, strides, paddings):
    bs = []

    for t in pooling_types:
        if t == "max":
            b = MaxPooling2D(window_lens, strides, paddings)
        elif t == 'avg':
            b = AvgPooling2D(window_lens, strides, paddings)
        else:
            raise ValueError

        bs.append(b)

    return Or(bs)

def MaybeSwap_fn(b1, b2):
    """Builds a module that has a parameter to swapping the order of modules
    passed as argument.
    """

    b = Or([
            Concat([b1, b2]),
            Concat([b2, b1])
        ])
    return b


text:
import os
import sys
import argparse
import numpy as np

from utils import *

def main(arguments):
    '''
    Read SICK official dataset and produce txt files for sent1, sent2 and similarity score
    '''
    parser = argparse.ArgumentParser(
        description=__doc__,
        formatter_class=argparse.ArgumentDefaultsHelpFormatter)
    parser.add_argument('--data_folder', help="location of folder with the sick files")
    parser.add_argument('--out_folder', help="location of the output folder")

    args = parser.parse_args(arguments)

    for split in ["train", "dev", "test"]:
        src_out = open(os.path.join(args.out_folder, "src-"+split+".txt"), "w")
        targ_out = open(os.path.join(args.out_folder, "targ-"+split+".txt"), "w")
        label_out = open(os.path.join(args.out_folder, "label-"+split+".txt"), "w")

        for line in open(os.path.join(args.data_folder, "SICK_"+split+".txt"),"r"):
            d = line.split("\t")
            label = d[4].strip().lower()
            premise = d[1].strip().lower()
            hypothesis = d[2].strip().lower()

            for punct in punctuations:
                if punct in premise:
                    premise = premise.replace(punct,"")
            #premise = [w for w in premise.split() if w not in stopwords]
            #premise = " ".join(premise)
            #premise = " ".join(lemmatize(" ".join(premise)))

            for punct in punctuations:
                if punct in hypothesis:
                    hypothesis = hypothesis.replace(punct,"")
            #hypothesis = [w for w in hypothesis.split() if w not in stopwords]
            #hypothesis = " ".join(hypothesis)
            #hypothesis = " ".join(lemmatize(" ".join(hypothesis)))

            src_out.write(premise + "\n")
            targ_out.write(hypothesis + "\n")
            label_out.write(label + "\n")

        src_out.close()
        targ_out.close()
        label_out.close()

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


text:
"""
Script to scrape airfoil profiles from http://airfoiltools.com/ with properties
for Re=200,000 and Ncrit=9.
Assumptions:
 - profile definition start at TE/SS and end at TE/PS and are continuous (Selig)
 - Consistent runs of profile simulations through XFoil

Run:
 > python3 airfoil_scrape.py --points
"""

################################################################################
# %% IMPORT PACKAGES
################################################################################

import os
import argparse
import requests
from bs4 import BeautifulSoup
import re
import numpy as np
import ray
import more_itertools as mit

################################################################################
# %% COMMAND LINE ARGUMENT PARSING
################################################################################

parser = argparse.ArgumentParser()
parser.add_argument("POINTS", help="number of points around circumference")
parser.add_argument("CPUS", help="number of CPUs to use for ray")
parser.add_argument("--cosine", help="use cosine spacing to refine LE", action="store_true")
args = parser.parse_args()

################################################################################
# %% CONSTANTS
################################################################################

BASE_URL = "http://airfoiltools.com"
POINTS = int(args.POINTS)
CPUS = int(args.CPUS)

################################################################################
# %% CREATE FOLDER IF NOT YET EXISTING
################################################################################

try:
    os.stat(f'01-prep-data/{POINTS}')
except:
    os.mkdir(f'01-prep-data/{POINTS}')

################################################################################
# %% INIT
################################################################################

ray.init(num_cpus=CPUS)

################################################################################
# %% SCAN FOR AIRFOILS
################################################################################

##### GRAB OVERVIEW PAGE
page = requests.get(BASE_URL+"/search/airfoils/index.html")
overview_soup = BeautifulSoup(page.text, 'html.parser')

##### FIND ALL LINKS
links = overview_soup.find_all("a", href=re.compile("airfoil/details"))
print(f'Found {len(links)} airfoil definitions')
links = [link.get('href') for link in links]

##### SPLIT INTO CPU-SIZED CHUNKS
chunk_links = [list(c) for c in mit.divide(CPUS, links)]

##### DEF REMOTE FUNCTION HERE!
@ray.remote
def scrape(links, BASE_URL):

    ##### INIT BATCH ARRAYS
    X_batch = np.zeros((0, POINTS, 2, 1), dtype=float)
    y_batch = np.zeros((0, 3), dtype=float)

    ##### LOOP OVER LINKS
    for link in links:

        ##### GRAB PROFILE PAGE
        page = requests.get(BASE_URL+link)
        airfoil_soup = BeautifulSoup(page.text, 'html.parser')

        ##### FIND DAT FILE CONTAINING PROFILE
        dat_link = airfoil_soup.find("a", href=re.compile("seligdat"))

        ##### GRAB PROFILE
        profile_page = requests.get(BASE_URL+dat_link.get('href'))
        profile_soup = BeautifulSoup(profile_page.text, 'html.parser')

        ##### PARSE PROFILE
        text = profile_soup.text.split('\n')
        profile_name = text[0]
        raw_profile = np.loadtxt(text[1:])

        ##### INTERPOLATE TO POINTS
        # RAW PROFILE
        fp1 = raw_profile[:, 0]
        fp2 = raw_profile[:, 1]

        # SCALE TO LOCAL DISTANCE
        delta = np.sqrt(np.diff(fp1)**2 + np.diff(fp2)**2)
        xp = delta.cumsum()
        xp = np.insert(xp, 0, 0)/xp[-1]

        # TARGET DISTRIBUTION
        if not args.cosine:
            ##### LINEAR
            x = np.linspace(0, 1, POINTS)
        elif args.cosine:
            ##### COSINE
            x = np.linspace(0,1, POINTS) + 0.1*np.sin(np.linspace(0, 2*np.pi, POINTS))

        # INTERPOLATE
        base_profile = np.array([np.interp(x, xp, fp1), np.interp(x, xp, fp2)]).T

        ##### SCALE BACK TO 0-1
        xmin = base_profile[:, 0].min()
        xmax = base_profile[:, 0].max()
        base_profile[:, 0] = (base_profile[:, 0]-xmin)/(xmax-xmin)
        base_profile[:, 1] = (base_profile[:, 1])/(xmax-xmin)

        ##### CALCULATE AREA
        area = 0.0
        for i in range(len(base_profile)-1):
            area += 0.5*(base_profile[i, 0]-base_profile[i+1, 0])*(base_profile[i, 1]+base_profile[i+1, 1])

        ##### FIND CL/CD DATA LINK
        data_link = airfoil_soup.find_all("a", href=re.compile("^/polar/details.*200000$"))

        # FOLLOW IF RE 200k DATA EXISTS AND AREA REASONABLE
        if data_link and area <= 0.2:
            data_page = requests.get(BASE_URL+data_link[0].get('href'))
            data_soup = BeautifulSoup(data_page.text, 'html.parser')

            ##### PARSE DATA TABLE
            data = []
            table = data_soup.find('table', attrs={'class':'tabdata'})
            rows = table.find_all('tr')
            for row in rows[1:]:
                cols = row.find_all('td')
                cols = [ele.text.strip() for ele in cols]
                data.append([ele for ele in cols if ele])
            data = np.array(data).astype(float)[:, 0:3]

            ##### INIT SAMPLE ARRAY
            X = np.zeros((len(data), POINTS, 2, 1), dtype=float)
            y = np.zeros((len(data), 3), dtype=float)

            ##### RUN THROUGH ANGLES AND BUILD PROFILES
            for sample in range(len(data)):

                ##### ANLGE IN RADIANS
                alpha = -data[sample, 0]*np.pi/180.0

                shift_profile = base_profile.copy()
                shift_profile[:, 0] = shift_profile[:, 0] - 0.5

                ##### ROTATE!
                profile_x = shift_profile[:, 0]*np.cos(alpha) - shift_profile[:, 1]*np.sin(alpha)
                profile_y = shift_profile[:, 0]*np.sin(alpha) + shift_profile[:, 1]*np.cos(alpha)

                ##### COMBINE
                profile_x = profile_x + 0.5
                profile = np.concatenate((profile_x.reshape(-1, 1), profile_y.reshape(-1, 1)), axis=1)

                ##### SET SAMPLE
                X[sample, :, :, 0] = profile
                y[sample, 0:2] = data[sample, 1:3]
                y[sample, 2] = area

            ##### COMBINE
            X_batch = np.concatenate((X_batch, X), axis=0)
            y_batch = np.concatenate((y_batch, y), axis=0)

    ##### RETURN!
    return X_batch, y_batch

################################################################################
# %% RUN IN PARALLEL!
################################################################################

object_id = [scrape.remote(chunk_links[batch], BASE_URL) for batch in range(CPUS)]
out = ray.get(object_id)

################################################################################
# %% SHUFFLE AND SAVE FILES
################################################################################

for i, batch in enumerate(out):

    ##### GET LOCAL BATCH
    X_batch = batch[0]
    y_batch = batch[1]

    ##### SHUFFLE
    idx = np.random.permutation(len(X_batch))
    X_batch = X_batch[idx]
    y_batch = y_batch[idx]

    ##### SAVE
    np.save(f'01-prep-data/{POINTS}/X_{i:03d}.npy', X_batch)
    np.save(f'01-prep-data/{POINTS}/y_{i:03d}.npy', y_batch)

max_stars_repo_name: markusmeingast/Airfoil-GAN

max_issues_repo_name: markusmeingast/Airfoil-GAN

max_forks_repo_name: markusmeingast/Airfoil-GAN


text:
import numpy as np
from scipy import optimize
import matplotlib.pyplot as plt
import math

# autoreload modules when code is run
#load_ext autoreload
#autoreload 2

# Question 1

# Construct a function which solves eq. (1)

# Following variables is known
m   = 1
v   = 10
e   = 0.3
t0  = 0.4
t1  = 0.1
k   = 0.4
w   = 0.5

# Now the model equations are defined

def x(m,w,l,t0,t1,k):
    """x equals the optimal consumption"""

    return m + w*l - (t0*w*l + t1*max(w*l-k,0))

def utility(l,c,e,v):
    return math.log(c) - v * ((l** (1+(1/e)))/(1+ (1/e)))

def labor_supply(l,e,v,m,w,t0,t1,k):
    c = x(m,w,l,t0,t1,k)
    return -utility(l,c,e,v)

def supply_problem(e,v,m,w,t0,t1,k):
    solution =  optimize.minimize_scalar(labor_supply,method='bounded', \
                bounds=(0,1), args=(e,v,m,w,t0,t1,k))
    l = solution.x
    c = x(m,w,l,t0,t1,k)
    u = utility(l,c,e,v)
    return l,c,u

l_star = supply_problem(e,v,m,w,t0,t1,k)[0]
c_star = supply_problem(e,v,m,w,t0,t1,k)[1]
u_star = supply_problem(e,v,m,w,t0,t1,k)[2]

print(f'Optimized labour supply is: {l_star:.2f}')
print(f'Optimized consumption is:  {c_star:.2f}')
print(f'Optimal utility is: {u_star:.2f}')

# Question 2

# A vector for w is created
N       = 10000
w_vec   = np.random.uniform(0.5, 1.5, size=N)

# empty lists for l and c is created
c_vec = np.empty(N)
l_vec = np.empty(N)

# a loop is generated to create optimal value of consumption \
# and labor supply given the wage

for i, w in enumerate(w_vec):
    opt = supply_problem(e,v,m,w_vec[i],t0,t1,k)
    c_vec[i] = opt[0]
    l_vec[i] = opt[1]

plt.style.use('grayscale')
plt.plot(w_vec,l_vec)
plt.plot(w_vec,c_vec)
plt.grid(True)
plt.show()

# Question 3

# A function to return the total tax revenue is created
np.random.seed(4600)
PopulationWageLarge = np.random.uniform(0.5,1.5,size=10000)

def TotalTax(PopulationWageVec,e,v,t0,t1,k,m):

#Return an array for individual tax payments
    N=len(PopulationWageVec)
    PopulationTaxes=np.zeros((N))

    for i,w in enumerate(PopulationWageVec):
        #Return optimal labour supply given optimize functions in Q1
        Ind_optimum=supply_problem(e,v,m,w,t0,t1,k)
        IndLabour=Ind_optimum[0]
        #Optimal invidual taxpayment with optimal labor
        PopulationTaxes[i]=t0*w*IndLabour+t1*max(w*IndLabour-k,0)
    #Sum
    TotTax=sum(PopulationTaxes)
    return TotTax

#Total tax functions with random uniform dis.
TotTax0 = TotalTax(PopulationWageLarge,e,v,t0,t1,k,m)
print(f'The total tax revenue is {TotTax0:.1f}')

# Question 4
# Frisch Elasticity changed from 0.3 to 0.1
e_new = 0.1

TotTax_e_new = TotalTax(PopulationWageLarge,e_new,v,t0,t1,k,m)
print(f'The total tax revenue is {TotTax_e_new:.1f}')

# Question 5
PopulationWageSmall = np.random.uniform(0.5,1.5,size=100)
PopulationWageMedium = np.random.uniform(0.5,1.5,size=1000)

def value_of_choice_tax(taxes,PopulationWage,e,v,m):

    t0 = taxes[0]
    t1 = taxes[1]
    k = taxes[2]
    return -TotalTax(PopulationWage,e,v,t0,t1,k,m)

#Finds the tax maximising values of the vector 'taxes'.
def taxOptimiser(PopulationWage,e,v,m):

    # 5.2.3 Call solver
    initial_guess = [0.8,0.6,0.5]
    sol_case3 = optimize.minimize(
        value_of_choice_tax,initial_guess,method='Nelder-Mead',
        args=(PopulationWage,e,v,m))

    t0Star=sol_case3.x[0]
    t1Star=sol_case3.x[1]
    kStar=sol_case3.x[2]

    #Print the solution
    print(f'Optimal standard income tax rate is {t0Star:.3f}')
    print(f'Optimal top bracet tax rate is {t1Star:.3f}')
    print(f'Optimal cut-off income is {kStar:.3f}')

    return[t0Star,t1Star,kStar]

print('Optimal taxe rates and revenue, with N=100')
[t0Star,t1Star,kStar]= taxOptimiser(PopulationWageSmall,e,v,m)
TotTaxSmall = TotalTax(PopulationWageLarge,e,v,t0Star,t1Star,kStar,m)
print(f'Total tax revenue = {TotTaxSmall:.2f}')

print('Optimal taxe rates and revenue, with N=1000')
[t0Star,t1Star,kStar]=taxOptimiser(PopulationWageMedium,e,v,m)
TotTaxMedium = TotalTax(PopulationWageLarge,e,v,t0Star,t1Star,kStar,m)
print(f'Total tax revenue = {TotTaxMedium:.2f}')

#print('Optimal taxe rates and revenue, with N=10000')
#[t0Star,t1Star,kStar]=taxOptimiser(PopulationWageLarge,e,v,m)
#TotTaxLarge = TotalTax(PopulationWageLarge,e,v,t0Star,t1Star,kStar,m)
#print(f'Total tax revenue = {TotTaxLarge:.2f}')

reason: import numpy,from scipy


text:
import numpy as np
from matplotlib import pyplot as plt
import animatplot as amp

def main():
    fig, ax = plt.subplots()
    ax.set_title("sin curve animation")
    ax.set_xlabel("x")
    ax.set_ylabel("y")
    xs = np.linspace(0, 1, 50)
    ts = np.linspace(0, 1, 20)
    Xs, Ts = np.meshgrid(xs, ts)
    Ys = np.sin(2 * np.pi * (Xs + Ts))
    block = amp.blocks.Line(Xs, Ys, axis=ax)
    timeline = amp.Timeline(ts, units='s', fps=20)
    anim = amp.Animation([block], timeline)
    anim.controls()
    anim.save_gif("sin_curve_with_contol")
    plt.show()

if __name__ == '__main__':
    main()

max_stars_repo_name: terasakisatoshi/pythonCodes

max_issues_repo_name: terasakisatoshi/pythonCodes

max_forks_repo_name: terasakisatoshi/pythonCodes


text:
"""Authors: Cody Baker and Ben Dichter."""
import pandas as pd
import numpy as np
from scipy.io import loadmat
from pathlib import Path
from lxml import etree as et
from datetime import datetime

from nwb_conversion_tools import NWBConverter
from nwb_conversion_tools import (
    NeuroscopeRecordingInterface,
    NeuroscopeSortingInterface,
)

from .yutalfpdatainterface import YutaLFPInterface, get_reference_elec
from .yutapositiondatainterface import YutaPositionInterface
from .yutabehaviordatainterface import YutaBehaviorInterface
from ..utils.neuroscope import get_clusters_single_shank, read_spike_clustering

def get_UnitFeatureCell_features(fpath_base, session_id, session_path, nshanks):
    """Load features from matlab file. Handle occasional mismatches."""
    cols_to_get = ("fineCellType", "region", "unitID", "unitIDshank", "shank")
    matin = loadmat(str(fpath_base / "DGProject/DG_all_6__UnitFeatureSummary_add.mat"), struct_as_record=False,)[
        "UnitFeatureCell"
    ][0][0]

    all_ids = []
    all_shanks = []
    for shankn in range(1, nshanks + 1):
        ids = np.unique(read_spike_clustering(session_path, shankn))
        ids = ids[~np.isin(ids, (0, 1))]
        all_ids.append(ids)
        all_shanks.append(np.ones(len(ids), dtype=int) * shankn)
    np.hstack(all_ids)
    np.hstack(all_shanks)
    clu_df = pd.DataFrame({"unitIDshank": np.hstack(all_ids), "shank": np.hstack(all_shanks)})

    this_file = matin.fname == session_id
    mat_df = pd.DataFrame({col: getattr(matin, col)[this_file].ravel() for col in cols_to_get})

    return pd.merge(clu_df, mat_df, how="left", on=("unitIDshank", "shank"))

class YutaNWBConverter(NWBConverter):
    """Primary conversion class for the SenzaiY dataset."""

    data_interface_classes = dict(
        NeuroscopeRecording=NeuroscopeRecordingInterface,
        NeuroscopeSorting=NeuroscopeSortingInterface,
        YutaLFP=YutaLFPInterface,
        YutaPosition=YutaPositionInterface,
        YutaBehavior=YutaBehaviorInterface,
    )

    def get_metadata(self):
        session_path = Path(self.data_interface_objects["NeuroscopeSorting"].source_data["folder_path"])
        subject_path = session_path.parent
        session_id = session_path.stem
        mouse_number = session_id[-9:-7]

        subject_xls = subject_path / f"DGProject/YM{mouse_number} exp_sheet.xlsx"
        hilus_csv_path = subject_path / "DGProject/early_session_hilus_chans.csv"
        if "-" in session_id:
            b = False
        else:
            b = True

        session_start = datetime.strptime(session_id[-6:], "%y%m%d")

        if subject_xls.is_file():
            subject_df = pd.read_excel(subject_xls)
            subject_data = dict()
            for key in [
                "genotype",
                "DOB",
                "implantation",
                "Probe",
                "Surgery",
                "virus injection",
                "mouseID",
            ]:
                names = subject_df.iloc[:, 0]
                if key in names.values:
                    subject_data[key] = subject_df.iloc[np.argmax(names == key), 1]
            if isinstance(subject_data["DOB"], datetime):
                age = str(session_start - subject_data["DOB"])
            else:
                age = None
        else:
            age = "unknown"
            subject_data = dict()
            print(f"Warning: no subject file detected for session {session_path}!")

        xml_filepath = session_path / f"{session_id}.xml"
        root = et.parse(str(xml_filepath)).getroot()

        shank_channels = [
            [int(channel.text) for channel in group.find("channels")]
            for group in root.find("spikeDetection").find("channelGroups").findall("group")
        ]

        all_shank_channels = np.concatenate(shank_channels)
        all_shank_channels.sort()
        nshanks = len(shank_channels)
        lfp_channel = get_reference_elec(subject_xls, hilus_csv_path, session_start, session_id, b=b)

        celltype_dict = {
            0: "unknown",
            1: "granule cells (DG) or pyramidal cells (CA3)  (need to use region info. see below.)",
            2: "mossy cell",
            3: "narrow waveform cell",
            4: "optogenetically tagged SST cell",
            5: "wide waveform cell (narrower, exclude opto tagged SST cell)",
            6: "wide waveform cell (wider)",
            8: "positive waveform unit (non-bursty)",
            9: "positive waveform unit (bursty)",
            10: "positive negative waveform unit",
        }

        df_unit_features = get_UnitFeatureCell_features(subject_path, session_id, session_path, nshanks)

        # there are occasional mismatches between the matlab struct
        # and the neuroscope files regions: 3: 'CA3', 4: 'DG'
        celltype_names = []
        for celltype_id, region_id in zip(df_unit_features["fineCellType"].values, df_unit_features["region"].values):
            if celltype_id == 1:
                if region_id == 3:
                    celltype_names.append("pyramidal cell")
                elif region_id == 4:
                    celltype_names.append("granule cell")
                else:
                    raise Exception("unknown type")
            elif not np.isfinite(celltype_id):
                celltype_names.append("missing")
            else:
                celltype_names.append(celltype_dict[celltype_id])

        sorting_electrode_groups = []
        for shankn in range(len(shank_channels)):
            df = get_clusters_single_shank(session_path, shankn + 1)
            for _, _ in df.groupby("id"):
                sorting_electrode_groups.append(f"shank{str(shankn+1)}")

        metadata = super().get_metadata()
        metadata["NWBFile"].update(
            session_start_time=session_start.astimezone(),
            session_id=session_id,
            institution="NYU",
            lab="Buzsaki",
        )
        metadata.update(Subject=dict(subject_id=session_id[:-7], species="Mus musculus", age=age))
        if "genotype" in subject_data:
            metadata["Subject"].update(genotype=subject_data["genotype"])
        metadata["Ecephys"]["Electrodes"].append(
            dict(
                name="theta_reference",
                description="This electrode was used to calculate LFP canonical bands.",
                data=list(all_shank_channels == lfp_channel),
            )
        )
        metadata.update(
            UnitProperties=[
                dict(
                    name="cell_type",
                    description="Name of cell type.",
                    data=celltype_names,
                ),
                dict(
                    name="global_id",
                    description="Global id for cell for entire experiment.",
                    data=df_unit_features["unitID"].values,
                ),
                dict(
                    name="shank_id",
                    description="0-indexed id of cluster of shank.",
                    # - 2 b/c the get_UnitFeatureCell_features removes 0 and 1 IDs from each shank
                    data=[x - 2 for x in df_unit_features["unitIDshank"].values],
                ),
                dict(
                    name="electrode_group",
                    description="The electrode group that each spike unit came from.",
                    data=sorting_electrode_groups,
                ),
            ]
        )

        return metadata

max_stars_repo_name: catalystneuro/buzsaki-lab-to-nwb

max_issues_repo_name: catalystneuro/buzsaki-lab-to-nwb


text:
import pytest
import numpy as np

import dowhy.api
from dowhy.causal_estimators.propensity_score_weighting_estimator import PropensityScoreWeightingEstimator
from .base import TestEstimator

from sklearn.linear_model import LinearRegression

class TestPropensityScoreMatchingEstimator(object):
    @pytest.mark.parametrize(["error_tolerance", "Estimator"],
                             [(0.05, PropensityScoreWeightingEstimator),
                              (0.1, PropensityScoreWeightingEstimator)])
    def test_average_treatment_effect(self, error_tolerance, Estimator):
        estimator_tester = TestEstimator(error_tolerance, Estimator)
        estimator_tester.average_treatment_effect_test()

    def test_pandas_api_discrete_cause_continuous_confounder(self):
        N = 10000
        data = dowhy.datasets.linear_dataset(beta=10,
                                             num_common_causes=1,
                                             num_instruments=1,
                                             num_samples=N,
                                             treatment_is_binary=False)
        X0 = np.random.normal(size=N)
        v = (np.random.normal(size=N) + X0).astype(int)
        y = data['ate']*v + X0 + np.random.normal()
        data['df']['v'] = v
        data['df']['X0'] = X0
        data['df']['y'] = y
        df = data['df'].copy()

        variable_types = {'v': 'd', 'X0': 'c', 'y': 'c'}
        outcome = 'y'
        cause = 'v'
        common_causes = 'X0'
        method = 'weighting'
        causal_df = df.causal.do(x=cause,
                                 variable_types=variable_types,
                                 outcome=outcome,
                                 method=method,
                                 common_causes=common_causes,
                                 proceed_when_unidentifiable=True)

        ate = (causal_df[causal_df.v == 1].mean() \
              - causal_df[causal_df.v == 0].mean())['y']
        error = np.abs(ate - data['ate'])
        error_tolerance = 0.05
        res = True if (error < data['ate'] * error_tolerance) else False
        print("Error in ATE estimate = {0} with tolerance {1}%. Estimated={2},True={3}".format(
            error, error_tolerance * 100, ate, data['ate'])
        )
        assert res

    def test_pandas_api_discrete_cause_discrete_confounder(self):
        N = 10000
        data = dowhy.datasets.linear_dataset(beta=10,
                                             num_common_causes=1,
                                             num_instruments=1,
                                             num_samples=N,
                                             treatment_is_binary=False)
        X0 = np.random.normal(size=N).astype(int)
        v = (np.random.normal(size=N) + X0).astype(int)
        y = data['ate'] * v + X0 + np.random.normal()
        data['df']['v'] = v
        data['df']['X0'] = X0
        data['df']['y'] = y
        df = data['df'].copy()

        variable_types = {'v': 'd', 'X0': 'd', 'y': 'c'}
        outcome = 'y'
        cause = 'v'
        common_causes = 'X0'
        method = 'weighting'
        causal_df = df.causal.do(x=cause,
                                 variable_types=variable_types,
                                 outcome=outcome,
                                 method=method,
                                 common_causes=common_causes,
                                 proceed_when_unidentifiable=True)

        ate = (causal_df[causal_df.v == 1].mean() \
              - causal_df[causal_df.v == 0].mean())['y']
        print('ate', ate)
        error = np.abs(ate - data['ate'])
        error_tolerance = 0.05
        res = True if (error < data['ate'] * error_tolerance) else False
        print("Error in ATE estimate = {0} with tolerance {1}%. Estimated={2},True={3}".format(
            error, error_tolerance * 100, ate, data['ate'])
        )
        assert res

    def test_pandas_api_continuous_cause_discrete_confounder(self):
        N = 1000
        data = dowhy.datasets.linear_dataset(beta=10,
                                             num_common_causes=1,
                                             num_instruments=1,
                                             num_samples=N,
                                             treatment_is_binary=False)
        X0 = np.random.normal(size=N).astype(int)
        v = np.random.normal(size=N) + X0
        y = data['ate'] * v + X0 + np.random.normal()
        data['df']['v'] = v
        data['df']['X0'] = X0
        data['df']['y'] = y
        df = data['df'].copy()

        variable_types = {'v': 'c', 'X0': 'd', 'y': 'c'}
        outcome = 'y'
        cause = 'v'
        common_causes = 'X0'
        method = 'weighting'
        causal_df = df.causal.do(x=cause,
                                 variable_types=variable_types,
                                 outcome=outcome,
                                 method=method,
                                 common_causes=common_causes,
                                 proceed_when_unidentifiable=True)

        ate = LinearRegression().fit(causal_df[['v']], causal_df['y']).coef_[0]
        print('ate', ate)
        error = np.abs(ate - data['ate'])
        error_tolerance = 0.05
        res = True if (error < data['ate'] * error_tolerance) else False
        print("Error in ATE estimate = {0} with tolerance {1}%. Estimated={2},True={3}".format(
            error, error_tolerance * 100, ate, data['ate'])
        )
        assert res

    def test_pandas_api_continuous_cause_continuous_confounder(self):
        N = 1000
        data = dowhy.datasets.linear_dataset(beta=10,
                                             num_common_causes=1,
                                             num_instruments=1,
                                             num_samples=N,
                                             treatment_is_binary=False)
        X0 = np.random.normal(size=N)
        v = np.random.normal(size=N) + X0
        y = data['ate'] * v + X0 + np.random.normal()
        data['df']['v'] = v
        data['df']['X0'] = X0
        data['df']['y'] = y
        df = data['df'].copy()

        variable_types = {'v': 'c', 'X0': 'c', 'y': 'c'}
        outcome = 'y'
        cause = 'v'
        common_causes = 'X0'
        method = 'weighting'
        causal_df = df.causal.do(x=cause,
                                 variable_types=variable_types,
                                 outcome=outcome,
                                 method=method,
                                 common_causes=common_causes,
                                 proceed_when_unidentifiable=True)

        ate = LinearRegression().fit(causal_df[['v']], causal_df['y']).coef_[0]
        print('ate', ate)
        error = np.abs(ate - data['ate'])
        error_tolerance = 0.05
        res = True if (error < data['ate'] * error_tolerance) else False
        print("Error in ATE estimate = {0} with tolerance {1}%. Estimated={2},True={3}".format(
            error, error_tolerance * 100, ate, data['ate'])
        )
        assert res

max_stars_repo_name: TechCityVentures/dowhy

max_issues_repo_name: TechCityVentures/dowhy

max_forks_repo_name: TechCityVentures/dowhy


text:
"""External Scorer for Beam Search Decoder."""
from __future__ import absolute_import
from __future__ import division
from __future__ import print_function

import os
import kenlm
import numpy as np

class Scorer(object):
    """External scorer to evaluate a prefix or whole sentence in
       beam search decoding, including the score from n-gram language
       model and word count.

    :param alpha: Parameter associated with language model. Don't use
                  language model when alpha = 0.
    :type alpha: float
    :param beta: Parameter associated with word count. Don't use word
                count when beta = 0.
    :type beta: float
    :model_path: Path to load language model.
    :type model_path: basestring
    """

    def __init__(self, alpha, beta, model_path, vocabulary):
        self._alpha = alpha
        self._beta = beta
        if not os.path.isfile(model_path):
            raise IOError("Invaid language model path: %s" % model_path)
        self._language_model = kenlm.LanguageModel(model_path)

    # n-gram language model scoring
    def _language_model_score(self, sentence):
        #log10 prob of last word
        log_cond_prob = list(
            self._language_model.full_scores(sentence, eos=False))[-1][0]
        return np.power(10, log_cond_prob)

    # word insertion term
    def _word_count(self, sentence):
        words = sentence.strip().split(' ')
        return len(words)

    # reset alpha and beta
    def reset_params(self, alpha, beta):
        self._alpha = alpha
        self._beta = beta

    # execute evaluation
    def __call__(self, sentence, log=False):
        """Evaluation function, gathering all the different scores
        and return the final one.

        :param sentence: The input sentence for evalutation
        :type sentence: basestring
        :param log: Whether return the score in log representation.
        :type log: bool
        :return: Evaluation score, in the decimal or log.
        :rtype: float
        """
        lm = self._language_model_score(sentence)
        word_cnt = self._word_count(sentence)
        if log == False:
            score = np.power(lm, self._alpha) * np.power(word_cnt, self._beta)
        else:
            score = self._alpha * np.log(lm) + self._beta * np.log(word_cnt)
        return score


text:
# main classification script

import keras
#import matplotlib.image as mpimg
import numpy as np
import itertools
import tensorflow as tf

from gan_dir_it import *

from imgaug import augmenters as iaa

from custom_resnet import ResNet18
from keras import backend as K
from keras.callbacks import TensorBoard, ModelCheckpoint, ReduceLROnPlateau
from keras.layers import Dense, Dropout
from keras.models import Model
from keras.optimizers import Adam
from keras.preprocessing import image

from math import ceil
from my_opts import opt
from os import listdir, makedirs, path, walk

if opt.mem_perc < 1.0:
  config = tf.ConfigProto()
  config.gpu_options.allow_growth = True
  config.gpu_options.per_process_gpu_memory_fraction = opt.mem_perc
  K.tensorflow_backend.set_session(tf.Session(config = config))

################################

# PART ONE: set names, directories etc

################################

data_dir = opt.data_dir
target_size = opt.target_size
stats_dir = data_dir + 'stats_train/'

name = 'bs_{}'.format(opt.batch_size)

dropout    = opt.dropout_prob > 0.0
fully_conn = opt.neurons > 0
augment    = opt.aug_prob > 0.0
fake_ratio = opt.fake_ratio

if (fake_ratio <= 0.0) or (fake_ratio > 2.0):
  print('Wrong fake ratio given. Should be greater than 0.0 and no greater than 2.0. Exiting.')
  import sys
  sys.exit()

if fake_ratio < 1.0:
  fake_ratio_str = '0' + str(int(opt.fake_ratio*100))
else:
  fake_ratio_str = str(int(opt.fake_ratio*100))

stats_path = stats_dir + 'stats_{}.txt'.format(fake_ratio_str)

if dropout:
  name += '_dr_{}'.format(int(opt.dropout_prob * 100))
if fully_conn:
  name += '_n_{}'.format(opt.neurons)

arch_name = name

name += '/{}'.format(fake_ratio_str)

if augment:
  name = 'aug_no_cont/aug_prob_{}/'.format(int(opt.aug_prob * 100)) + name
else:
  name = 'no_aug/' + name
name += '/'

makedirs(name)
print(name)

if opt.weights_dir is None:
  weights_path = name
  print('yes')
else:
  print('no')
  weights_path = name + opt.weights_dir
  if not opt.weights_dir.endswith('/'):
    weights_path += '/'
  makedirs(weights_path)

best_weights_path = weights_path + 'best_weights_{}.h5'.format(arch_name)
print(best_weights_path)

if opt.log_dir is not None:
  log_dir = name + opt.log_dir
  makedirs(log_dir)

batch_size = opt.batch_size

################################

# PART TWO: preprocess input data, set pipelines etc

################################

def data_mean_std():
    names = []
    for root, dirs, files in walk(data_dir):
        if "train" in root:
            for file in files:
                names.append(root + '/' + file)

    array = []
    for name in names:
        img = image.load_img(name, grayscale = True, target_size = target_size)

        img_array = np.asarray(img)
        array.append(img_array)

    array = np.asarray(array)
    return array.mean(), array.std()

aug_prob = opt.aug_prob
seq = iaa.Sequential([
        iaa.Fliplr(aug_prob),
	iaa.Sometimes(aug_prob, iaa.Multiply((0.9, 1.1))),
        iaa.Sometimes(aug_prob, iaa.Affine(
                  scale={"x": (0.9, 1.1), "y": (0.9, 1.1)},  # scale images to 80-120% of their size
                  translate_percent={"x": (-0.05, 0.05), "y": (-0.05, 0.05)},  # translate by -5 to +5 percent (per axis)
                  rotate=(-5, 5),  # rotate by -5 to +5 degrees
                  cval=0,
                  mode='constant'
                  ))
      ])

#mean, std = data_mean_std()
#mean, std = 35.31741725575326, 44.92010876082367
stats_file = open(stats_path, 'r')
stats_lst = stats_file.read().split('\n')
mean = float(stats_lst[0])
std = float(stats_lst[1])

print(mean, std)
def normalize(x):
    return (x - mean) / std

def norm_and_aug(x):
    return seq.augment_image(normalize(x))

train_path = path.join(data_dir, "train_{}".format(fake_ratio_str))
valid_path = path.join(data_dir, "valid")

#gen = image.ImageDataGenerator(featurewise_center = True, featurewise_std_normalization = True)
if aug_prob > 0.0:
  train_gen = image.ImageDataGenerator(preprocessing_function = norm_and_aug)
#  train_gen = GANImageDataGenerator(preprocessing_function = norm_and_aug)
else:
  train_gen = image.ImageDataGenerator(preprocessing_function = normalize)
#  train_gen = GANImageDataGenerator(preprocessing_function = normalize)

valid_gen = image.ImageDataGenerator(preprocessing_function = normalize)

train_batches = train_gen.flow_from_directory(train_path, target_size = target_size, color_mode = 'grayscale', batch_size = batch_size)
valid_batches = valid_gen.flow_from_directory(valid_path, target_size = target_size, color_mode = 'grayscale', batch_size = batch_size)

################################

# PART THREE: load model

################################

tb = TensorBoard(log_dir = log_dir, batch_size = batch_size, write_graph = True)

lr_reducer = ReduceLROnPlateau(monitor='loss', factor=np.sqrt(0.1), patience=10, cooldown=0, min_lr=1e-5)

model_chkp = ModelCheckpoint(best_weights_path,
                            # monitor='val_quadratic_kappa_score',
                            monitor='val_acc',
                            save_best_only=True, save_weights_only=True, mode='max')

resnet = ResNet18()
#resnet = ResNet50()
#model = keras.applications.resnet50.ResNet50(include_top=True, weights='imagenet', input_tensor=None, input_shape=None, pooling=None, classes=1000)
#model.summary()

#resnet.layers.pop()
#for layer in resnet.layers:
#    layer.trainable = False

x = resnet.layers[-1].output
if dropout:
  x = Dropout(opt.dropout_prob)(x)
if fully_conn:
  x = Dense(opt.neurons, activation = 'relu')(x)
  if dropout:
    x = Dropout(opt.dropout_prob)(x)

x = Dense(2, activation = 'softmax')(x)

ft_resnet = Model(resnet.input, x)
ft_resnet.compile(optimizer = Adam(lr = 0.001), loss = 'categorical_crossentropy', metrics = ['accuracy'])

print("model loaded")

################################

# PART FOUR: train model

################################

epochs = opt.epochs
train_samples = len(listdir(train_path + '/ad/')) + len(listdir(train_path + '/nor/'))
valid_samples = len(listdir(valid_path + '/ad/')) + len(listdir(valid_path + '/nor/'))

print(train_samples, valid_samples)

train_steps = ceil(train_samples/batch_size)
valid_steps = ceil(valid_samples/batch_size)

print (train_steps, valid_steps)

ft_resnet.fit_generator(train_batches, steps_per_epoch=train_steps,
                        epochs=epochs, callbacks = [tb, model_chkp, lr_reducer],
                        validation_data=valid_batches, validation_steps = valid_steps)

continue_model_path = weights_path + '{}_{}_eps.h5'.format(arch_name, epochs)
ft_resnet.save(continue_model_path)

max_issues_repo_name: nikkkkhil/GAN-MRI-data-augmentation

max_forks_repo_name: nikkkkhil/GAN-MRI-data-augmentation