def test_energy():
    # make sure that energy as computed by ssvm is the same as by lp
    np.random.seed(0)
    for inference_method in ["lp", "ad3"]:
        found_fractional = False
        crf = EdgeFeatureGraphCRF(n_states=3,
                                  inference_method=inference_method,
                                  n_edge_features=2)
        while not found_fractional:
            x = np.random.normal(size=(7, 8, 3))
            edge_list = make_grid_edges(x, 4, return_lists=True)
            edges = np.vstack(edge_list)
            edge_features = edge_list_to_features(edge_list)
            x = (x.reshape(-1, 3), edges, edge_features)

            unary_params = np.random.normal(size=(3, 3))
            pw1 = np.random.normal(size=(3, 3))
            pw2 = np.random.normal(size=(3, 3))
            w = np.hstack([unary_params.ravel(), pw1.ravel(), pw2.ravel()])
            res, energy = crf.inference(x, w, relaxed=True, return_energy=True)
            found_fractional = np.any(np.max(res[0], axis=-1) != 1)

            psi = crf.psi(x, res)
            energy_svm = np.dot(psi, w)

            assert_almost_equal(energy, -energy_svm)
            if not found_fractional:
                # exact discrete labels, test non-relaxed version
                res, energy = crf.inference(x, w, relaxed=False,
                                            return_energy=True)
                psi = crf.psi(x, res)
                energy_svm = np.dot(psi, w)

                assert_almost_equal(energy, -energy_svm)

def test_query_resolution():
    bbox = (-180, -90, 180, 90)
    init_res = (4.5, 6.7)
    query = mapnik.Query(mapnik.Box2d(*bbox), init_res)
    r = query.resolution
    assert_almost_equal(r[0], init_res[0], places=7)
    assert_almost_equal(r[1], init_res[1], places=7)

def test_keep_top_k_epitopes():
    arg_parser = make_variant_sequences_arg_parser()
    args = arg_parser.parse_args([
        "--vcf", data_path("b16.f10/b16.f10.Phip.vcf"),
        "--bam", data_path("b16.f10/b16.combined.sorted.bam"),
    ])
    reads_generator = allele_reads_generator_from_args(args)
    variants = variant_collection_from_args(args)
    keep_k_epitopes = 3
    core_logic = VaxrankCoreLogic(
        reads_generator=reads_generator,
        mhc_predictor=random_binding_predictor,
        variants=variants,
        vaccine_peptide_length=15,
        padding_around_mutation=5,
        min_alt_rna_reads=1,
        min_variant_sequence_coverage=1,
        variant_sequence_assembly=True,
        max_vaccine_peptides_per_variant=1,
        num_mutant_epitopes_to_keep=keep_k_epitopes)
    ranked_list = core_logic.ranked_vaccine_peptides()

    for variant, vaccine_peptides in ranked_list:
        vaccine_peptide = vaccine_peptides[0]
        eq_(keep_k_epitopes, len(vaccine_peptide.mutant_epitope_predictions))
        # recompute the expected score, make sure the top-k argument from ranked_vaccine_peptides()
        # propagated as expected
        mutant_epitope_score = sum(
            p.logistic_epitope_score() for p in vaccine_peptide.mutant_epitope_predictions)
        assert_almost_equal(mutant_epitope_score, vaccine_peptide.mutant_epitope_score)

def test_avoid_nearby_birds():
    regression_data=yaml.load(open(os.path.join(os.path.dirname(__file__),'avoid_nearby_boids.yml')))
    boid_data=regression_data["before"]
    avoid_nearby_boids(boid_data)
    for after,before in zip(regression_data["after"],boid_data):
        for after_value,before_value in zip(after,before): 
            assert_almost_equal(after_value,before_value,delta=0.01)

def test_bad_boids_regression():
    regression_data=yaml.load(open(os.path.join(os.path.dirname(__file__),'fixture.yml')))
    boid_data=regression_data["before"]
    update_boids(boid_data)
    for after,before in zip(regression_data["after"],boid_data):
        for after_value,before_value in zip(after,before): 
            assert_almost_equal(after_value,before_value,delta=0.01)

def test_logging():
    iris = load_iris()
    X, y = iris.data, iris.target

    X_ = [(np.atleast_2d(x), np.empty((0, 2), dtype=np.int)) for x in X]
    Y = y.reshape(-1, 1)

    X_train, X_test, y_train, y_test = train_test_split(X_, Y, random_state=1)
    _, file_name = mkstemp()

    pbl = GraphCRF(n_features=4, n_states=3, inference_method=inference_method)
    logger = SaveLogger(file_name)
    svm = NSlackSSVM(pbl, C=100, n_jobs=1, logger=logger)
    svm.fit(X_train, y_train)

    score_current = svm.score(X_test, y_test)
    score_auto_saved = logger.load().score(X_test, y_test)

    alt_file_name = file_name + "alt"
    logger.save(svm, alt_file_name)
    logger.file_name = alt_file_name
    logger.load()
    score_manual_saved = logger.load().score(X_test, y_test)

    assert_less(.97, score_current)
    assert_less(.97, score_auto_saved)
    assert_less(.97, score_manual_saved)
    assert_almost_equal(score_auto_saved, score_manual_saved)

def _check_marginal_samples_match_scores(server, row, fi):
    row = loom.query.protobuf_to_data_row(row.diff)
    row[fi] = None
    to_sample = [i == fi for i in range(len(row))]
    samples = server.sample(to_sample, row, SAMPLE_COUNT)
    val = samples[0][fi]
    base_score = server.score(row)
    if isinstance(val, bool) or isinstance(val, int):
        probs_dict = {}
        samples = [sample[fi] for sample in samples]
        for sample in set(samples):
            row[fi] = sample
            probs_dict[sample] = numpy.exp(
                server.score(row) - base_score)
        if len(probs_dict) == 1:
            assert_almost_equal(probs_dict[sample], 1., places=SCORE_PLACES)
            return
        if min(probs_dict.values()) < MIN_CATEGORICAL_PROB:
            return
        gof = discrete_goodness_of_fit(samples, probs_dict, plot=True)
    elif isinstance(val, float):
        probs = numpy.exp([
            server.score(sample) - base_score
            for sample in samples
        ])
        samples = [sample[fi] for sample in samples]
        gof = density_goodness_of_fit(samples, probs, plot=True)
    assert_greater(gof, MIN_GOODNESS_OF_FIT)

    def test_remove_momentum_rescale_energy_openmm(self):
        # don't actually need to do everything with OpenMM, but do need to
        # add units
        u_vel = old_div(u.nanometer, u.picosecond)
        u_mass = old_div(u.dalton, u.AVOGADRO_CONSTANT_NA)
        u_energy = old_div(u.kilojoule_per_mole, u.AVOGADRO_CONSTANT_NA)

        velocities = \
                np.array([[1.5, -1.0], [-1.0, 2.0], [0.25, -1.0]]) * u_vel
        masses = np.array([1.0, 1.5, 4.0]) * u_mass
        new_vel = self.openmm_modifier._remove_linear_momentum(
            velocities=velocities,
            masses=masses
        )
        new_momenta = new_vel * masses[:, np.newaxis]
        total_momenta = sum(new_momenta, new_momenta[0])
        assert_array_almost_equal(total_momenta,
                                  np.array([0.0]*2) * u_vel * u_mass)

        new_vel = self.openmm_modifier._rescale_kinetic_energy(
            velocities=velocities,
            masses=masses,
            double_KE=20.0 * u_energy
        )
        new_momenta = new_vel * masses[:, np.newaxis]
        total_momenta = sum(new_momenta, new_momenta[0])
        zero_energy = 0.0 * u_energy
        new_ke = sum(sum(new_momenta * new_vel, zero_energy), zero_energy)
        # tests require that the linear momentum be 0, and KE be correct
        assert_array_almost_equal(total_momenta,
                                  np.array([0.0]*2) * u_vel * u_mass)
        assert_equal(new_ke.unit, (20.0 * u_energy).unit)
        assert_almost_equal(new_ke._value, (20.0 * u_energy)._value)

    def test_call(self):
        new_toy_snap = self.toy_modifier(self.toy_snapshot)
        assert_array_almost_equal(new_toy_snap.coordinates,
                                  self.toy_snapshot.coordinates)
        new_vel = new_toy_snap.velocities
        old_vel = self.toy_snapshot.velocities
        same_vel = [np.allclose(new_vel[i], old_vel[i]) 
                    for i in range(len(new_vel))]
        assert_equal(Counter(same_vel), Counter({True: 2, False: 1}))
        for new_v, old_v in zip(new_vel, old_vel):
            assert_almost_equal(sum([v**2 for v in new_v]),
                                sum([v**2 for v in old_v]))

        new_omm_snap = self.openmm_modifier(self.openmm_snap)
        n_atoms = len(self.openmm_snap.coordinates)
        assert_array_almost_equal(new_omm_snap.coordinates,
                                  self.openmm_snap.coordinates)
        new_vel = new_omm_snap.velocities
        old_vel = self.openmm_snap.velocities
        same_vel = [np.allclose(new_vel[i], old_vel[i]) 
                    for i in range(len(new_vel))]
        same_vel = [np.allclose(new_vel[i], old_vel[i]) 
                    for i in range(len(new_vel))]
        assert_equal(Counter(same_vel), Counter({True: n_atoms-1, False: 1}))
        u_vel_sq = (old_div(u.nanometers, u.picoseconds))**2
        for new_v, old_v in zip(new_vel, old_vel):
            assert_almost_equal(
                sum([(v**2).value_in_unit(u_vel_sq) for v in new_v]),
                sum([(v**2).value_in_unit(u_vel_sq) for v in old_v])
            )

def test_initial_value():
    m = Minuit(func3, pedantic=False, x=1., y=2., error_x=3., print_level=0)
    assert_almost_equal(m.args[0], 1.)
    assert_almost_equal(m.args[1], 2.)
    assert_almost_equal(m.values['x'], 1.)
    assert_almost_equal(m.values['y'], 2.)
    assert_almost_equal(m.errors['x'], 3.)

 def test_intercepts_from_means_with_true_normalization(self):
     expected_intercepts = self.true_intercepts_series
     true_normalization = ['m3', 2.0]
     calc_intercepts, calc_mean = wf.intercepts_from_means(
         self.data, true_normalization, self.true_loadings_series)
     assert_series_equal(calc_intercepts, expected_intercepts)
     assert_almost_equal(calc_mean, 0.0, places=2)

def check_consistency_in_mesh1d(macro_grid, element_orders, var_list):
    my_mesh1d = Mesh1D(macro_grid, element_orders, var_list)

    # Testing mesh attributes
    numpy.testing.assert_allclose(my_mesh1d.macro_grid, macro_grid)
    numpy.testing.assert_array_equal(my_mesh1d.element_orders, element_orders)
    assert_equal(my_mesh1d.variables, var_list)
    assert_equal(my_mesh1d.dof, (sum(element_orders) + 1) * len(var_list))

    # Testing list of elements
    for idx_var, var in enumerate(var_list):
        for idx_el, element in enumerate(my_mesh1d.elem):
            if idx_var == 0:
                mesh_pos_1 = my_mesh1d.gm[idx_el][element.pos[var]][0]
                mesh_pos_2 = sum(my_mesh1d.elem[el].order for el in range(idx_el))
                numpy.testing.assert_array_equal(mesh_pos_1, mesh_pos_2)
            else:
                mesh_pos_1 = my_mesh1d.gm[idx_el][element.pos[previous_var]] + my_mesh1d.dof_1v
                mesh_pos_2 = my_mesh1d.gm[idx_el][element.pos[var]]
                numpy.testing.assert_array_equal(mesh_pos_1, mesh_pos_2)
        previous_var = var

    # Testing linear and quadratic integration
    domain_integral_a = sum(el.w_1v.dot(el.x_1v) for el in my_mesh1d.elem)
    domain_integral_b = sum(el.w_1v.dot(el.x_1v ** 2) for el in my_mesh1d.elem)
    assert_almost_equal(domain_integral_a, (my_mesh1d.macro_grid[-1]**2 - my_mesh1d.macro_grid[0]**2) / 2)
    assert_almost_equal(domain_integral_b, (my_mesh1d.macro_grid[-1]**3 - my_mesh1d.macro_grid[0]**3) / 3)

 def test_fit_background_C(self):
     self.m["B_K"].active = False
     self.m.fit_background()
     nt.assert_almost_equal(self.m["Offset"].offset.value,
                            1.71212121212)
     nt.assert_false(self.m["B_K"].active)
     nt.assert_true(self.m["C_K"].active)

def test_generate_og_receptive_field():
    xpixels = 500 # simulated screen width
    ypixels = 500 # simulated screen height
    ppd = 1 # simulated visual angle
    scale_factor = 1.0 # simulated stimulus resampling rate
    distance = 5 # standard deviations to compute gauss out to
    xcenter = 0 # x coordinate of the pRF center
    ycenter = 0 # y coordinate of the pRF center
    sigma = 1 # width of the pRF
    
    test_value = 6 # this is the sum of a gaussian given 1 ppd 
                   # and a 1 sigma prf centered on (0,0)
                   
    # generate the visuotopic coordinates
    dx,dy = generate_coordinate_matrices(xpixels,
                                         ypixels,
                                         ppd,
                                         scale_factor)
    
    # generate a pRF at (0,0) and 1 sigma wide
    rf = generate_og_receptive_field(xcenter, ycenter, sigma, dx, dy)
    
    # divide by integral
    rf /= 2 * np.pi * sigma ** 2
    
    # compare the volume of the pRF to a known value
    nt.assert_almost_equal(np.sum(rf),1)

 def test_x(self):
     s = session.query(Spot).get(1)
     assert_almost_equal(float(session.scalar(s.spot_location.x)), -88.594586159235689)
     s = session.query(Spot).filter(and_(Spot.spot_location.x < 0, Spot.spot_location.y > 42)).all()
     ok_(s is not None)
     assert_almost_equal(float(session.scalar(functions.x(WKTSpatialElement('POINT(-88.3655256496815 43.1402866687898)', geometry_type=Point.name)))),
                         -88.3655256496815)

    def test_call_with_linear_momentum_fix(self):
        toy_modifier = SingleAtomVelocityDirectionModifier(
            delta_v=[1.0, 2.0],
            subset_mask=[1, 2],
            remove_linear_momentum=True
        )
        new_toy_snap = toy_modifier(self.toy_snapshot)
        velocities = new_toy_snap.velocities
        momenta = velocities * new_toy_snap.masses[:, np.newaxis]
        assert_array_almost_equal(sum(momenta), np.array([0.0]*2))
        double_ke = sum(sum(momenta * velocities))
        assert_almost_equal(double_ke, 86.0)

        u_vel = old_div(u.nanometer, u.picosecond)
        u_mass = old_div(u.dalton, u.AVOGADRO_CONSTANT_NA)

        openmm_modifier = SingleAtomVelocityDirectionModifier(
            delta_v=1.2*u_vel,
            remove_linear_momentum=False
        )
        new_openmm_snap = openmm_modifier(self.openmm_snap)
        velocities = new_openmm_snap.velocities
        momenta = velocities * new_openmm_snap.masses[:, np.newaxis]
        zero_momentum = 0 * u_vel * u_mass
        total_momenta = sum(momenta, zero_momentum)
        assert_array_almost_equal(total_momenta,
                                  np.array([0.0]*3) * u_vel * u_mass)

 def test_y(self):
     s = session.query(Spot).get(1)
     assert_almost_equal(float(session.scalar(s.spot_location.y)), 42.9480095987261)
     s = session.query(Spot).filter(and_(Spot.spot_location.y < 0, Spot.spot_location.y > 42)).all()
     ok_(s is not None)
     assert_almost_equal(float(session.scalar(functions.y(WKTSpatialElement('POINT(-88.3655256496815 43.1402866687898)', geometry_type=Point.name)))),
                         43.1402866687898)

def test_mean_frequency1():
    """basic: Test mean_frequency 1"""

    import efel
    import numpy

    stim_start = 500.0
    stim_end = 900.0

    data = numpy.loadtxt('testdata/basic/mean_frequency_1.txt')

    time = data[:, 0]
    voltage = data[:, 1]

    trace = {}

    trace['T'] = time
    trace['V'] = voltage
    trace['stim_start'] = [stim_start]
    trace['stim_end'] = [stim_end]

    features = ['mean_frequency']

    feature_values = \
        efel.getFeatureValues(
            [trace],
            features)
    nt.assert_almost_equal(feature_values[0]['mean_frequency'], 15.2858453)

 def test_length(self):
     l = session.query(Lake).get(1)
     r = session.query(Road).get(1)
     s = session.query(Spot).get(1)
     assert_almost_equal(session.scalar(l.lake_geom.length), 0.30157858985653774)
     assert_almost_equal(session.scalar(r.road_geom.length), 0.8551694164147895)
     ok_(not session.scalar(s.spot_location.length))

def test_ap_amplitude_from_voltagebase1():
    """basic: Test AP_amplitude_from_voltagebase 1"""

    import efel
    import numpy

    stim_start = 500.0
    stim_end = 900.0

    data = numpy.loadtxt('testdata/basic/mean_frequency_1.txt')

    time = data[:, 0]
    voltage = data[:, 1]

    trace = {}

    trace['T'] = time
    trace['V'] = voltage
    trace['stim_start'] = [stim_start]
    trace['stim_end'] = [stim_end]

    features = ['AP_amplitude_from_voltagebase',
                'peak_voltage', 'voltage_base']

    feature_values = \
        efel.getFeatureValues(
            [trace],
            features)

    voltage_base = feature_values[0]['voltage_base'][0]
    for peak_voltage, ap_amplitude_from_voltagebase in zip(
            feature_values[0]['peak_voltage'],
            feature_values[0]['AP_amplitude_from_voltagebase']):
        nt.assert_almost_equal(peak_voltage - voltage_base,
                               ap_amplitude_from_voltagebase)