def test():
    dummy_data = {
        'PRES': ma.masked_array([1.0, 100, 200, 300, 500, 5000]),
        'TEMP': ma.masked_array([27.44, 14.55, 11.96, 11.02, 7.65, 2.12]),
        'PSAL': ma.masked_array([35.71, 35.50, 35.13, 35.02, 34.72, 35.03])
        }
    features = {
            'tukey53H': ma.masked_array([0, 0, 0.3525000000000009,
                0.35249999999999915, 0, 0],
                mask=[True, True, False, False, True, True]),
            'tukey53H_norm': ma.masked_array([0, 0, 0.07388721803621254,
                0.07388721803621218, 0, 0],
                mask = [True,  True, False, False, True, True])
            }
    flags = {'tukey53H_norm': np.array([0, 0, 1, 1, 0, 0], dtype='i1')}

    cfg = {
            'l': 5,
            'threshold': 6,
            'flag_good': 1,
            'flag_bad': 4
            }

    y = Tukey53H(dummy_data, 'TEMP', cfg)
    y.test()

    assert type(y.features) is dict
    for f in y.features:
        assert ma.allclose(y.features[f], features[f])
    for f in y.flags:
        assert ma.allclose(y.flags[f], flags[f])

 def test_fully_qualified_DDD_SSS_variable(self):
     # this should use the test variable a_test_SSS_variable_DDD_SSS
     expr = "opus_core.tests.a_test_squid_variable_42_clam"
     storage = StorageFactory().get_storage('dict_storage')
     storage.write_table(
         table_name='tests',
         table_data={
             "a_dependent_variable":array([1,5,10]),
             "id":array([1,3,4])
             }
         )
     dataset = Dataset(in_storage=storage, in_table_name='tests', id_name="id", dataset_name="tests")
     result = dataset.compute_variables([expr])
     should_be = array([10,50,100])
     self.assert_(ma.allclose(result, should_be, rtol=1e-6), "Error in test_fully_qualified_DDD_SSS_variable")
     # check that the access methods for the variable all return the correct values
     name = VariableName(expr)
     self.assertEqual(name.get_package_name(), 'opus_core', msg="bad value for package")
     self.assertEqual(name.get_dataset_name(), 'tests', msg="bad value for dataset")
     self.assertEqual(name.get_short_name(), 'a_test_squid_variable_42_clam', msg="bad value for shortname")
     self.assertEqual(name.get_alias(), 'a_test_squid_variable_42_clam', msg="bad value for alias")
     self.assertEqual(name.get_autogen_class(), None, msg="bad value for autogen_class")
     # test that the variable can now also be accessed using its short name in an expression
     result2 = dataset.compute_variables(['a_test_squid_variable_42_clam'])
     self.assert_(ma.allclose(result2, should_be, rtol=1e-6), "Error in accessing a_test_squid_variable_42_clam")

 def test_fully_qualified_variable(self):
     # this tests an expression consisting of a fully-qualified variable
     expr = "opus_core.test_agent.income_times_2"
     storage = StorageFactory().get_storage("dict_storage")
     storage.write_table(table_name="test_agents", table_data={"income": array([1, 5, 10]), "id": array([1, 3, 4])})
     dataset = Dataset(in_storage=storage, in_table_name="test_agents", id_name="id", dataset_name="test_agent")
     result = dataset.compute_variables([expr])
     should_be = array([2, 10, 20])
     self.assert_(ma.allclose(result, should_be, rtol=1e-6), "Error in test_fully_qualified_variable")
     # check that expr is in the cache of known expressions
     # (normally we shouldn't be accessing this private field, but just this once ...)
     cache = VariableName._cache
     self.assert_(expr in cache, msg="did not find expr in cache")
     # check that the access methods for the variable all return the correct values
     name = VariableName(expr)
     self.assertEqual(name.get_package_name(), "opus_core", msg="bad value for package")
     self.assertEqual(name.get_dataset_name(), "test_agent", msg="bad value for dataset")
     self.assertEqual(name.get_short_name(), "income_times_2", msg="bad value for shortname")
     self.assertEqual(name.get_alias(), "income_times_2", msg="bad value for alias")
     self.assertEqual(name.get_autogen_class(), None, msg="bad value for autogen_class")
     # test that the variable can now also be accessed using its short name in an expression
     result2 = dataset.compute_variables(["income_times_2"])
     self.assert_(ma.allclose(result2, should_be, rtol=1e-6), "Error in accessing a_test_variable")
     # check that the cache uses the variable name with whitespace removed
     oldsize = len(cache)
     expr_with_spaces = "opus_core . test_agent. income_times_2  "
     name2 = VariableName(expr_with_spaces)
     newsize = len(cache)
     self.assertEqual(oldsize, newsize, msg="caching error")
     self.assert_(expr_with_spaces not in cache, msg="caching error")
     self.assertEqual(expr_with_spaces, name2.get_expression(), msg="caching error")
     self.assertEqual(name2.get_short_name(), "income_times_2", msg="bad value for shortname")

def isEqual(left, right, eps=None, masked_equal=True):
  ''' This function checks if two numpy arrays or scalars are equal within machine precision, and returns a scalar logical. '''
  diff_type = "Both arguments to function 'isEqual' must be of the same class!"
  if isinstance(left,np.ndarray):
    # ndarray
    if not isinstance(right,np.ndarray): raise TypeError(diff_type)
    if not left.dtype==right.dtype:
      right = right.astype(left.dtype) # casting='same_kind' doesn't work...
    if np.issubdtype(left.dtype, np.inexact): # also catch float32 etc
      if eps is None: return ma.allclose(left, right, masked_equal=masked_equal)
      else: return ma.allclose(left, right, masked_equal=masked_equal, atol=eps)
    elif np.issubdtype(left.dtype, np.integer) or np.issubdtype(left.dtype, np.bool):
      return np.all( left == right ) # need to use numpy's all()
  elif isinstance(left,(float,np.inexact)):
    # numbers
    if not isinstance(right,(float,np.inexact)): raise TypeError(diff_type)
    if eps is None: eps = 100.*floateps # default
    if ( isinstance(right,float) or isinstance(right,float) ) or left.dtype.itemsize == right.dtype.itemsize: 
      return np.absolute(left-right) <= eps
    else:
      if left.dtype.itemsize < right.dtype.itemsize: right = left.dtype.type(right)
      else: left = right.dtype.type(left)
      return np.absolute(left-right) <= eps  
  elif isinstance(left,(int,bool,np.integer,np.bool)):
    # logicals
    if not isinstance(right,(int,bool,np.integer,np.bool)): raise TypeError(diff_type)
    return left == right
  else: raise TypeError(left)

    def test_deletion_of_jobs_with_one_characteristics(self):
        dataset_pool = DatasetPool(storage=self.storage, package_order=["washtenaw","urbansim", "opus_core"])
        gridcell_set = dataset_pool.get_dataset('gridcell')
        event_set = self._create_job_deletion_event_set_with_characteristics()
        jobs = dataset_pool.get_dataset("job")
        AgentEventModel().run(gridcell_set, event_set, jobs, 2000, dataset_pool)
        number_of_jobs_of_sector_1 = gridcell_set.compute_variables("urbansim.gridcell.number_of_jobs_of_sector_1", 
                                                                    dataset_pool=dataset_pool)
        number_of_jobs_of_sector_2 = gridcell_set.compute_variables("urbansim.gridcell.number_of_jobs_of_sector_2", 
                                                                    dataset_pool=dataset_pool)
        number_of_jobs_of_sector_4 = gridcell_set.compute_variables("urbansim.gridcell.number_of_jobs_of_sector_4", 
                                                                    dataset_pool=dataset_pool)
        # the model should remove 2 jobs of sector 1 from gridcell 1, 
        #                         5 jobs of sector 1 from gridcell 5
        self.assert_(ma.allclose(number_of_jobs_of_sector_1, array( [2,4,4,4,0,4,4,4,4,4]))) 
        # other sectors don't change
        self.assert_(ma.allclose(number_of_jobs_of_sector_2, array( 10 * [3]))) 
        self.assert_(ma.allclose(number_of_jobs_of_sector_4, array( 10 * [3]))) 

        AgentEventModel().run(gridcell_set, event_set, jobs, 2001, dataset_pool)
        number_of_jobs_of_sector_1 = gridcell_set.compute_variables("urbansim.gridcell.number_of_jobs_of_sector_1", 
                                                                    dataset_pool=dataset_pool)
        number_of_jobs_of_sector_2 = gridcell_set.compute_variables("urbansim.gridcell.number_of_jobs_of_sector_2", 
                                                                    dataset_pool=dataset_pool)
        number_of_jobs_of_sector_4 = gridcell_set.compute_variables("urbansim.gridcell.number_of_jobs_of_sector_4", 
                                                                    dataset_pool=dataset_pool)
        # the model should remove 2 jobs of sector 2 from gridcell 5, 
        #                         1 job of sector 1 from gridcell 1,
        #                       all jobs of sector 2 from gridcell 2
        #                        70% jobs of sector 2 from gridcell 3
        self.assert_(ma.allclose(number_of_jobs_of_sector_1, array( [1,4,4,4,0,4,4,4,4,4]))) 
        self.assert_(ma.allclose(number_of_jobs_of_sector_2, array( [3, 0, 1, 3, 1, 3, 3, 3, 3, 3]))) 

    def test_land_price_model(self):
        storage = StorageFactory().get_storage('dict_storage')

        gridcell_set_table_name = 'gridcell_set'
        storage.write_table(
            table_name=gridcell_set_table_name,
            table_data={
                "percent_residential_within_walking_distance":array([30, 0, 90, 100]),
                "gridcell_year_built":array([2002, 1968, 1880, 1921]),
                "fraction_residential_land":array([0.5, 0.1, 0.3, 0.9]),
                "residential_land_value":array([0, 0, 0, 0]),
                "nonresidential_land_value":array([0, 0, 0, 0]),
                "development_type_id":array(  [1, 1,  1, 1]),
                "grid_id": array([1,2,3,4])
                }
            )
        gridcell_set = GridcellDataset(in_storage=storage, in_table_name=gridcell_set_table_name)

        specification = EquationSpecification(variables=(
            "percent_residential_within_walking_distance",
            "gridcell_year_built", "constant"),
            coefficients=("PRWWD", "YB", "constant"))
        coefficients = Coefficients(names=("constant", "PRWWD", "YB"), values=(10.0, -0.0025, 0.0001))
        lp = LandPriceModel(filter=None, debuglevel=3)
        lp.run(specification, coefficients, gridcell_set)
        result1 = gridcell_set.get_attribute("residential_land_value")
        result2 = gridcell_set.get_attribute("nonresidential_land_value")
        self.assertEqual(ma.allclose(result1, array([12482.124,  2681.723,  6367.914, 18708.617]), rtol=1e-3), True)
        self.assertEqual(ma.allclose(result2, array([12482.124,  24135.510, 14858.466,  2078.735]), rtol=1e-3), True)

    def test_same_distribution_after_household_subtraction(self):
        """Using the control_totals and no marginal characteristics,
        subtract households and ensure that the distribution within each group stays the same
        """
        annual_household_control_totals_data = {
            "year": array([2000]),
            "total_number_of_households": array([20000])
            }

        storage = StorageFactory().get_storage('dict_storage')

        storage.write_table(table_name='hh_set', table_data=self.households_data)
        hh_set = HouseholdDataset(in_storage=storage, in_table_name='hh_set')

        storage.write_table(table_name='hct_set', table_data=annual_household_control_totals_data)
        hct_set = ControlTotalDataset(in_storage=storage, in_table_name='hct_set', what="household", id_name="year")

        storage.write_table(table_name='hc_set', table_data=self.household_characteristics_for_ht_data)
        hc_set = HouseholdCharacteristicDataset(in_storage=storage, in_table_name='hc_set')

        model = HouseholdTransitionModel()
        model.run(year=2000, household_set=hh_set, control_totals=hct_set, characteristics=hc_set)

        #check that there are indeed 20000 total households after running the model
        results = hh_set.size()
        should_be = [20000]
        self.assertEqual(ma.allclose(should_be, results, rtol=1e-1),
                         True, "Error, should_be: %s, but result: %s" % (should_be, results))

        #check that the distribution of households in each group is the same as before running the model
        results = self.get_count_all_groups(hh_set)
        should_be = [6000.0/33000.0*20000.0, 2000.0/33000.0*20000.0, 3000.0/33000.0*20000.0, 4000.0/33000.0*20000.0,
                     2000.0/33000.0*20000.0, 5000.0/33000.0*20000.0, 3000.0/33000.0*20000.0, 8000.0/33000.0*20000.0]
        self.assertEqual(ma.allclose(results, should_be, rtol=0.05),
                         True, "Error, should_be: %s,\n but result: %s" % (should_be, results))

    def test_safely_divide_two_attributes(self):
        from opus_core.datasets.dataset_pool import DatasetPool
        
        storage = StorageFactory().get_storage('dict_storage')        
        storage.write_table(
            table_name='tests',
            table_data={
                'id': array([1,2,3,4]),
                'numerator': array([1,2,3,0]),
                'denominator': array([2.,0.,2.,0.]),
            }
        )
        
        dataset_pool = DatasetPool(package_order=['opus_core'],
                                   storage=storage)
        test = dataset_pool.get_dataset('test')
        variable = Variable()
        variable.set_dataset(test)

        result = variable.safely_divide_two_attributes('opus_core.test.numerator',
                                                       'opus_core.test.denominator')
        self.assert_(ma.allclose(array([.5, 0, 1.5, 0]), result))
        
        result = variable.safely_divide_two_attributes('opus_core.test.numerator',
                                                       'opus_core.test.denominator', 
                                                        value_for_divide_by_zero=-1.0)
        self.assert_(ma.allclose(array([.5, -1., 1.5, -1.]), result))

 def test_interaction_set_component(self):
     # test a fully-qualified variable that applies to a component of an interaction set
     expr = "opus_core.test_agent.income_times_2"
     storage = StorageFactory().get_storage('dict_storage')
     storage.write_table(
         table_name='test_agents', 
         table_data={'id': array([1, 2, 3]), 'income': array([1, 20, 500])}
         )
     storage.write_table(
         table_name='test_locations', 
         table_data={'id': array([1,2]), 'cost': array([1000, 2000])}
         )
     dataset_pool = DatasetPool(package_order=['opus_core'], storage=storage)
     test_agent_x_test_location = dataset_pool.get_dataset('test_agent_x_test_location')
     result = test_agent_x_test_location.compute_variables(expr, dataset_pool=dataset_pool)
     should_be = array([[2, 2], [40, 40], [1000, 1000]])
     self.assert_(ma.allclose(result, should_be, rtol=1e-6), msg = "Error in " + expr)
     # test that the interaction set now has this as an attribute
     result2 = test_agent_x_test_location.get_attribute('income_times_2')
     self.assert_(ma.allclose(result2, should_be, rtol=1e-6), msg = "Error in " + expr)
     # test that the variable can now also be accessed using its short name
     result3 = test_agent_x_test_location.compute_variables(['income_times_2'])
     self.assert_(ma.allclose(result3, should_be, rtol=1e-6), msg = "Error in " + expr)
     # even though we're using this with an interaction set, the dataset name for expr
     # should be the name of the component set (since that's the only one mentioned in expr)
     name = VariableName(expr)
     self.assertEqual(name.get_dataset_name(), 'test_agent', msg="bad value for dataset")

 def test_join_by_rows(self):
     storage = StorageFactory().get_storage('dict_storage')
     
     storage.write_table(
         table_name='dataset1', 
         table_data={    
             'id':array([2,4,6,8]), 
             'attr':array([4,7,2,1])
             }
         )
         
     storage.write_table(
         table_name='dataset2',
         table_data={
             'id':array([1,5,9]), 
             'attr':array([55,66,100])
             }
         )
     
     ds1 = Dataset(in_storage=storage, in_table_name='dataset1', id_name='id')
     ds2 = Dataset(in_storage=storage, in_table_name='dataset2', id_name='id')
     
     ds1.join_by_rows(ds2)
     self.assert_(ma.allclose(ds1.get_attribute('attr'), array([4,7,2,1,55,66,100])))
     self.assert_(ma.allclose(ds2.get_attribute('attr'), array([55,66,100])))

    def test_same_distribution_after_household_addition(self):
        """Using the control_totals and no marginal characteristics,
        add households and ensure that the distribution within each group stays the same
        """

        annual_household_control_totals_data = {
            "year": array([2000, 2000]),
            "total_number_of_households": array([20000, 30000]),
            "large_area_id": array([1,2])
            }

        storage = StorageFactory().get_storage('dict_storage')

        storage.write_table(table_name = 'hh_set', table_data = self.households_data)
        hh_set = HouseholdDataset(in_storage=storage, in_table_name='hh_set')
        
        storage.write_table(table_name = 'hct_set', table_data = annual_household_control_totals_data)
        hct_set = ControlTotalDataset(in_storage=storage, in_table_name='hct_set', what="household")

        storage.write_table(table_name = 'hc_set', table_data = self.household_characteristics_for_ht_data)
        hc_set = HouseholdCharacteristicDataset(in_storage=storage, in_table_name='hc_set')

        model = RegionalHouseholdTransitionModel()
        model.run(year=2000, household_set=hh_set, control_totals=hct_set, characteristics=hc_set)

        #check that there are 20000 (area 1) and 30000 (area 2) total households after running the model
        areas = hh_set.get_attribute("large_area_id")
        results = array([0,0])
        for iarea in [0,1]:
            results[iarea] = where(areas == [1,2][iarea])[0].size
        should_be = [20000, 30000]
        self.assertEqual(ma.allclose(should_be, results, rtol=1e-1),
                         True, "Error, should_be: %s, but result: %s" % (should_be, results))

        #check that the number of unplaced households is exactly the number of new households created
        results = where(hh_set.get_attribute("grid_id")<=0)[0].size
        should_be = [17000]
        self.assertEqual(ma.allclose(should_be, results, rtol=1e-1),
                         True, "Error, should_be: %s, but result: %s" % (should_be, results))

        #check that the distribution of households in each group and each area is the same as before running the model
        results = self.get_count_all_groups(hh_set)
        should_be = array([
                    # area 1 
                     3000.0/16500.0*20000.0, 1000.0/16500.0*20000.0, 1500.0/16500.0*20000.0, 2000.0/16500.0*20000.0,
                     1000.0/16500.0*20000.0, 2500.0/16500.0*20000.0, 1500.0/16500.0*20000.0, 4000.0/16500.0*20000.0,
                     # area 2
                     3000.0/16500.0*30000.0, 1000.0/16500.0*30000.0, 1500.0/16500.0*30000.0, 2000.0/16500.0*30000.0,
                     1000.0/16500.0*30000.0, 2500.0/16500.0*30000.0, 1500.0/16500.0*30000.0, 4000.0/16500.0*30000.0])
        self.assertEqual(ma.allclose(results, should_be, rtol=0.1),
                         True, "Error, should_be: %s, but result: %s" % (should_be, results))
        # check the types of the attributes
        self.assertEqual(hh_set.get_attribute("age_of_head").dtype, int32,
                         "Error in data type of the new household set. Should be: int32, is: %s" % str(hh_set.get_attribute("age_of_head").dtype))
        self.assertEqual(hh_set.get_attribute("income").dtype, int32,
                         "Error in data type of the new household set. Should be: int32, is: %s" % str(hh_set.get_attribute("income").dtype))
        self.assertEqual(hh_set.get_attribute("persons").dtype, int8,
                         "Error in data type of the new household set. Should be: int8, is: %s" % str(hh_set.get_attribute("persons").dtype))

 def test_simple_model_with_outcome_values(self):
     m = SimpleModel()
     m.run(self.dataset,  outcome_attribute='iniattr', outcome_values=zeros(10)-1)
     self.assertEqual(ma.allclose(self.dataset.get_attribute('iniattr'), array(10*[-1])), True)
     self.assertEqual('iniattr' in self.dataset.get_primary_attribute_names(), True)
     # run with filter
     m.run(self.dataset,  outcome_attribute='iniattr', outcome_values=arange(10)+1, dataset_filter='dataset.attribute>1000')
     expected = array([1, 2, -1, -1, -1, -1, 7, -1, -1, -1])
     self.assertEqual(ma.allclose(self.dataset.get_attribute('iniattr'), expected), True)

def test_estimate_anomaly():
    f1 = estimate_anomaly(dummy_features,
            {'spike': dummy_params['spike']})
    f2 = estimate_anomaly(pd.DataFrame(dummy_features),
            {'spike': dummy_params['spike']})

    assert ma.allclose(f1, f2)
    assert ma.allclose(f1,
            ma.masked_values([-999, 0.0, -5.797359001920061,
                -57.564627324851145, -999, -9.626760611162082], -999))

    def test_fave_scale(self) :
        hanning.hanning_smooth(self.Data)
        rebin_freq.rebin(self.Data, 2.)
        cal_scale.scale_by_cal(self.Data, False, True, False)
        data = self.Data.data

        self.assertTrue(ma.allclose(ma.mean(data[:,0,0,:] -
                                              data[:,0,1,:], -1), 1.0))
        self.assertTrue(ma.allclose(ma.mean(data[:,3,0,:] -
                                              data[:,3,1,:], -1), 1.0))

    def test_controlling_with_three_marginal_characteristics(self):
        """Controlling with all three possible marginal characteristics in this example, age_of_head, income, and persons,
        this would partition the 8 groups into the same 8 groups, and with a control total specified for each group, we must
        ensure that the control totals for each group exactly meet the specifications.
        """

        #IMPORTANT: marginal characteristics grouping indices have to start at 0!
        annual_household_control_totals_data = {
            "year": array(8*[2000]),
            #"age_of_head": array(4*[0] + 4*[1]),
            "age_of_head_min": array([ 0, 0, 0, 0, 50, 50, 50, 50]),
            "age_of_head_max": array([49,49,49,49,100,100,100,100]),
            #"income": array(2*[0] + 2*[1] + 2*[0] + 2*[1]),
            "income_min": array([    0,    0,40000,40000,    0,    0,40000,40000]),
            "income_max": array([39999,39999,   -1,   -1,39999,39999,   -1,   -1]),
            #"persons": array([0,1,0,1,0,1,0,1]),
            "persons_min": array([0, 3,0, 3,0, 3,0, 3]),
            "persons_max": array([2,-1,2,-1,2,-1,2,-1]),
            "total_number_of_households": array([4000, 5000, 1000, 3000, 0, 6000, 3000, 8000])
            }
        ##size of columns was not even, removed last element of min and max
        #household_characteristics_for_ht_data = {
            #"characteristic": array(2*['age_of_head'] + 2*['income'] + 2*['persons']),
            #"min": array([0, 50, 0, 40000, 0, 3]),
            #"max": array([49, 100, 39999, -1, 2, -1]) 
            #}
        storage = StorageFactory().get_storage('dict_storage')

        storage.write_table(table_name='hh_set', table_data=self.households_data)
        hh_set = HouseholdDataset(in_storage=storage, in_table_name='hh_set')

        storage.write_table(table_name='hct_set', table_data=annual_household_control_totals_data)
        hct_set = ControlTotalDataset(in_storage=storage, in_table_name='hct_set', what='household', id_name=[])

        #storage.write_table(table_name='hc_set', table_data=household_characteristics_for_ht_data)
        #hc_set = HouseholdCharacteristicDataset(in_storage=storage, in_table_name='hc_set')

        # unplace some households
        where10 = where(hh_set.get_attribute("grid_id")<>10)[0]
        hh_set.modify_attribute(name="grid_id", data=zeros(where10.size), index=where10)

        model = TransitionModel(hh_set, control_total_dataset=hct_set)
        model.run(year=2000, target_attribute_name="total_number_of_households", reset_dataset_attribute_value={'grid_id':-1})

        #check that there are indeed 33000 total households after running the model
        results = hh_set.size()
        should_be = [30000]
        self.assertEqual(ma.allclose(should_be, results, rtol=1e-1),
                         True, "Error, should_be: %s, but result: %s" % (should_be, results))

        #check that the number of households in each group exactly match the control totals specified
        results = self.get_count_all_groups(hh_set)
        should_be = [4000, 5000, 1000, 3000, 0, 6000, 3000, 8000]
        self.assertEqual(ma.allclose(results, should_be),
                         True, "Error, should_be: %s, but result: %s" % (should_be, results))

 def test_safely_divide_two_arrays(self):
     result = Variable().safely_divide_two_arrays(array([10,20,30,0]).astype(int8), array([2,0,2,0]).astype(int8))
     self.assert_(ma.allclose(array([5,0,15,0]), result))
     # Types are done correctly
     self.assertEqual(result.dtype.type, int8)
     
     result = Variable().safely_divide_two_arrays(array([1,2,3,0]), array([2.,0.,2.,0.]))
     self.assert_(ma.allclose(array([.5, 0, 1.5, 0]), result))
     
     result = Variable().safely_divide_two_arrays(array([1,2,3,0]), array([2.,0.,2.,0.]), value_for_divide_by_zero=-1.)
     self.assert_(ma.allclose(array([.5, -1., 1.5, -1.]), result))

    def test_controlling_with_one_marginal_characteristic(self):
        """Using the age_of_head as a marginal characteristic, which would partition the 8 groups into two larger groups
        (those with age_of_head < 40 and >= 40), ensure that the control totals are met and that the distribution within
        each large group is the same before and after running the model
        """

        #IMPORTANT: marginal characteristics grouping indices have to start at 0!
        #i.e. below, there is one marg. char. "age_of_head". here we indicate that the first "large group" (groups 1-4),
        #consisting of those groups with age_of_head < 40 should total 25000 households after running this model for one year,
        #and the second large group, those groups with age_of_head > 40, should total 15000 households
        annual_household_control_totals_data = {
            "year": array([2000, 2000]),
            "age_of_head": array([0,1]),
            "total_number_of_households": array([25000, 15000])
            }

        storage = StorageFactory().get_storage('dict_storage')

        storage.write_table(table_name='hh_set', table_data=self.households_data)
        hh_set = HouseholdDataset(in_storage=storage, in_table_name='hh_set')

        storage.write_table(table_name='hct_set', table_data=annual_household_control_totals_data)
        hct_set = ControlTotalDataset(in_storage=storage, in_table_name='hct_set', what='household', id_name=['year' ,'age_of_head'])

        storage.write_table(table_name='hc_set', table_data=self.household_characteristics_for_ht_data)
        hc_set = HouseholdCharacteristicDataset(in_storage=storage, in_table_name='hc_set')

        storage.write_table(table_name='prs_set', table_data=self.person_data)
        prs_set = PersonDataset(in_storage=storage, in_table_name='prs_set')
        
        model = HouseholdTransitionModel()
        model.run(year=2000, person_set=prs_set, household_set=hh_set, control_totals=hct_set, characteristics=hc_set)

        #check that there are indeed 40000 total households after running the model
        results = hh_set.size()
        should_be = [40000]
        self.assertEqual(ma.allclose(should_be, results, rtol=1e-1),
                         True, "Error, should_be: %s, but result: %s" % (should_be, results))

        #check that the total number of households within first four groups increased by 10000
        #and that the total number of households within last four groups decreased by 3000
        results = self.get_count_all_groups(hh_set)
        should_be = [25000, 15000]
        self.assertEqual(ma.allclose([sum(results[0:4]), sum(results[4:8])], should_be, rtol=1e-1),
                         True, "Error, should_be: %s, but result: %s" % (should_be, results))

        #check that the distribution of households within groups 1-4 and 5-8 are the same before and after
        #running the model, respectively

        should_be = [6000.0/15000.0*25000.0, 2000.0/15000.0*25000.0, 3000.0/15000.0*25000.0, 4000.0/15000.0*25000.0,
                     2000.0/18000.0*15000.0, 5000.0/18000.0*15000.0, 3000.0/18000.0*15000.0, 8000.0/18000.0*15000.0]
        self.assertEqual(ma.allclose(results, should_be, rtol=0.05),
                         True, "Error, should_be: %s, but result: %s" % (should_be, results))

    def test_same_distribution_after_household_addition(self):
        """Using the control_totals and no marginal characteristics,
        add households and ensure that the distribution within each group stays the same
        """

        annual_household_control_totals_data = {
            "year": array([2000]),
            "total_number_of_households": array([50000])
            }

        storage = StorageFactory().get_storage('dict_storage')

        storage.write_table(table_name='hh_set', table_data=self.households_data)
        hh_set = HouseholdDataset(in_storage=storage, in_table_name='hh_set')

        storage.write_table(table_name='prs_set', table_data=self.person_data)
        prs_set = PersonDataset(in_storage=storage, in_table_name='prs_set')
        
        storage.write_table(table_name='hct_set', table_data=annual_household_control_totals_data)
        hct_set = ControlTotalDataset(in_storage=storage, in_table_name='hct_set', what="household", id_name="year")

        storage.write_table(table_name='hc_set', table_data=self.household_characteristics_for_ht_data)
        hc_set = HouseholdCharacteristicDataset(in_storage=storage, in_table_name='hc_set')

        model = HouseholdTransitionModel()
        model.run(year=2000, person_set=prs_set, household_set=hh_set, control_totals=hct_set, characteristics=hc_set)

        #check that there are indeed 50000 total households after running the model
        results = hh_set.size()
        should_be = [50000]
        self.assertEqual(ma.allclose(should_be, results, rtol=1e-1),
                         True, "Error, should_be: %s, but result: %s" % (should_be, results))

        #check that the number of unplaced households is exactly the number of new households created
        results = where(hh_set.get_attribute("building_id")<=0)[0].size
        should_be = [17000]
        self.assertEqual(ma.allclose(should_be, results, rtol=1e-1),
                         True, "Error, should_be: %s, but result: %s" % (should_be, results))

        #check that the distribution of households in each group is the same as before running the model
        results = self.get_count_all_groups(hh_set)
        should_be = array([6000.0/33000.0*50000.0, 2000.0/33000.0*50000.0, 3000.0/33000.0*50000.0, 4000.0/33000.0*50000.0,
                     2000.0/33000.0*50000.0, 5000.0/33000.0*50000.0, 3000.0/33000.0*50000.0, 8000.0/33000.0*50000.0])
        self.assertEqual(ma.allclose(results, should_be, rtol=0.05),
                         True, "Error, should_be: %s, but result: %s" % (should_be, results))
        # check the types of the attributes
        self.assertEqual(hh_set.get_attribute("age_of_head").dtype, int32,
                         "Error in data type of the new household set. Should be: int32, is: %s" % str(hh_set.get_attribute("age_of_head").dtype))
        self.assertEqual(hh_set.get_attribute("income").dtype, int32,
                         "Error in data type of the new household set. Should be: int32, is: %s" % str(hh_set.get_attribute("income").dtype))
        self.assertEqual(hh_set.get_attribute("persons").dtype, int8,
                         "Error in data type of the new household set. Should be: int8, is: %s" % str(hh_set.get_attribute("persons").dtype))

    def test_addition_of_jobs(self):
        dataset_pool = DatasetPool(storage=self.storage, package_order=["washtenaw","urbansim", "opus_core"])
        gridcell_set = dataset_pool.get_dataset('gridcell')
        event_set = self._create_simple_job_addition_event_set()
        jobs = dataset_pool.get_dataset("job")
        AgentEventModel().run(gridcell_set, event_set, jobs, 2000, dataset_pool)
        number_of_jobs = gridcell_set.compute_variables("urbansim.gridcell.number_of_jobs", dataset_pool=dataset_pool)
        # the model should add 5 jobs to gridcell 5 and 20 jobs to gridcell 10
        self.assert_(ma.allclose(number_of_jobs, array( [10,10,10,10,15,10,10,10,10,30]))) 

        AgentEventModel().run(gridcell_set, event_set, jobs, 2001, dataset_pool)
        number_of_jobs = gridcell_set.compute_variables("urbansim.gridcell.number_of_jobs", dataset_pool=dataset_pool)
        # the model should add another 3 jobs to gridcell 5
        self.assert_(ma.allclose(number_of_jobs, array( [10,10,10,10,18,10,10,10,10,30])))

    def test_deletion_of_jobs(self):
        dataset_pool = DatasetPool(storage=self.storage, package_order=["washtenaw","urbansim", "opus_core"])
        gridcell_set = dataset_pool.get_dataset('gridcell')
        event_set = self._create_simple_deletion_event_set()

        DeletionEventModel().run(gridcell_set, event_set, 2000, dataset_pool)
        number_of_jobs = gridcell_set.compute_variables("urbansim.gridcell.number_of_jobs", dataset_pool=dataset_pool)
        # the model should remove 5 jobs from gridcell 5 and all jobs from gridcell 10
        self.assert_(ma.allclose(number_of_jobs, array( [10,10,10,10,5,10,10,10,10,0]))) 

        DeletionEventModel().run(gridcell_set, event_set, 2001, dataset_pool)
        number_of_jobs = gridcell_set.compute_variables("urbansim.gridcell.number_of_jobs", dataset_pool=dataset_pool)
        # the model should remove another 3 jobs from gridcell 5
        self.assert_(ma.allclose(number_of_jobs, array( [10,10,10,10,2,10,10,10,10,0])))