def ismember(element, array, rows=False):
    """Check if element is member of array"""

    if rows:
        return sp.any([sp.all(array[x, :] == element) for x in range(array.shape[0])])
    else:
        return sp.all([element[i] in array for i in element.shape[0]])

def plot_saved():
    import matplotlib
    import matplotlib.pyplot as plt

    data = sp.memmap(get_filename(), dtype=sp.float64, mode='r', shape=(N_samp + 2, N_steps / res_every + 1, N))
    
    exa = data[-1]
    
    data = data[:-2]
    
    fins = sp.array([d[plot_res] for d in data if not sp.all(d[plot_res] == 0)])
    nsam = len(fins)
    print("Samples:", nsam)
    av = fins.sum(axis=0) / nsam

    av_var = 1./(nsam - 1) / nsam * sp.sum((fins/nsam - av)**2, axis=0) 
    av_e1 = av + sp.sqrt(av_var)
    av_e2 = av - sp.sqrt(av_var)

    
    plt.figure()
    pav = plt.plot(av, 'k-')[0]
    plt.plot(av_e1, 'k--')
    plt.plot(av_e2, 'k--')

    if not sp.all(exa[-1] == 0):
        pexa = plt.plot(exa[-1], 'r-')[0]
        plt.legend([pexa, pav], ["Density matrix", "Sample average"])
    
    plt.ylim((-1, 1))
    
    plt.xlabel(r"$n$")
    plt.ylabel(r"$\langle \sigma^z_n \rangle$")
    
    plt.show()

 def test_block_mesh_dict(self, data_field_class):
     field = data_field_class()
     offsets = data_field_class()
     #
     params = {
         'convertToMeters': '0.000010000',
         'numbersOfCells': '(5 10 15)',
         'cellExpansionRatios': 'simpleGrading (1 2 3)',
         #
         'boundary.left.type': 'empty',
         'boundary.right.type': 'empty',
         'boundary.top.type': 'wall',
         'boundary.bottom.type': 'wall',
         'boundary.front.type': 'wall',
         'boundary.back.type': 'wall'
     }
     mesh = OpenFoam.BlockMeshDict(field, avg_fact=10.0, mesh_params=params, offset_field=offsets)
     mesh._edges = ['placeholder']
     mesh._mergePatchPairs = ['placeholder']
     mesh.write_foam_file(TEMP_DIR, overwrite=True)
     #
     # attempting to overwrite existing mesh file
     with pytest.raises(FileExistsError):
         mesh.write_mesh_file(TEMP_DIR, overwrite=False)
     #
     # writing out a symmetry plane
     mesh.write_symmetry_plane(TEMP_DIR, overwrite=True)
     #
     # testing generation of a thresholded mesh
     mesh.generate_threshold_mesh(min_value=9, max_value=90)
     assert sp.all(mesh.data_map[0, :] == 0)
     assert sp.all(mesh.data_map[9, :] == 0)
     assert len(mesh._blocks) == 80

def test_hessian():
    """ Test the the numerical hessian for the rosenbrock function 
    with and without explicit gradient and a function of a higher order
    """
    def F(x):
        return x[0]**3+x[1]**3+x[2]**3+x[0]**2 *x[1]**2 *x[2]**2
    def F_hessian(x):
        return array([ [6.*x[0] + 2.*x[1]**2*x[2]**2, 4.*x[0]*x[1]*x[2]**2, 4.*x[0]*x[1]**2 *x[2]] , 
                        [4.*x[0]*x[1]*x[2]**2, 6.*x[1] + 2.*x[0]**2*x[2]**2, 4.*x[0]**2 *x[1]*x[2]] ,
                        [4.*x[0]*x[1]**2 *x[2], 4.*x[0]**2 *x[1]*x[2],6.*x[2] + 2.*x[0]**2*x[1]**2 ]])
    opt1 = p.OptimizationProblem(rosen)
    opt2 = p.OptimizationProblem(rosen, rosen_der)
    opt3 = p.OptimizationProblem(F,3)
    for i in range(-3,3):
        for j in range(-3,3):
            x = array([i, j, 3], dtype=double)
            k  = opt1.hessian(x) - rosen_hess(x)
            kk  = opt2.hessian(x) - rosen_hess(x)
            kkk  = opt3.hessian(x) - F_hessian(x)
            print k, abs(k) <1e-2
            print kk, abs(kk) <1e-2
            print kkk, abs(kkk) <1e-2
            assert all( abs(k) <1e-2 )
            assert all( abs(kk) <1e-2 )
            assert all( abs(kkk) <1e-2 )

def _assert_files_equal(expected_path, actual_path):
    def o(f):
        if f.endswith('.gz'):
            return gzip.open(f, 'rb')
        elif f.endswith('.hdf5'):
            return h5py.File(f, 'r')
        elif f.endswith('.ps'):
            return open(f, 'r')
        elif f.endswith('.pickle'):
            return open(f, 'rb')
        else:
            return open(f, 'rb')

    with o(expected_path) as e:
        with o(actual_path) as a:
            if expected_path.endswith('.hdf5'):
                _compare_hdf5(e, a)
            elif expected_path.endswith('.ps'):
                _compare_ps(e, a)
            elif expected_path.endswith('.pickle'):
                ta = pickle.load(a, encoding='latin1')
                te = pickle.load(e, encoding='latin1')
                if len(ta) == 0 or isinstance(ta[0], sp.int64):
                    assert sp.all(ta == te)
                else:
                    assert sp.all([_compare_gene(_[0], _[1]) for _ in zip(ta, te)])
            elif os.path.basename(expected_path).startswith('test_results'):
                _assert_files_equal_testing(e, a)
            else:
                assert e.read() == a.read(), 'actual and expected content differ!\nactual path: %s\nexpected path: %s\n' % (expected_path, actual_path)

 def test_multiphase_init(self):
     m = op.phases.MultiPhase(network=self.net, phases=[self.air,
                                                        self.water])
     assert sp.all(m['pore.occupancy.all'] == 0.0)
     assert sp.all(m['throat.occupancy.all'] == 0.0)
     assert self.air.name in m.settings['phases']
     assert self.water.name in m.settings['phases']

def _calc_pcs(weight_dict, sids, nts, snps):
    num_nt_issues = 0
    num_snps_used = 0
    num_indivs = snps.shape[1]
    pcs = sp.zeros((num_indivs, 2))

    for snp_i, sid in enumerate(sids):
        try:
            d = weight_dict[sid]
        except:
            continue
        nt = nts[snp_i]
        snp = snps[snp_i]
        pc_weights = sp.array([d['pc1w'], d['pc2w']])
        pc_weights.shape = (1, 2)
        af = d['mean_g'] / 2.0
        if sp.all([nt[1], nt[0]] == d['nts']):
            # print 'Flip sign'
            pc_weights = -pc_weights
            af = 1 - af
        elif not sp.all(nt == d['nts']):
            num_nt_issues += 1
            continue
        mean_g = 2 * af
        sd_g = sp.sqrt(af * (1 - af))
        # "Normalizing" the SNPs with the given allele frequencies
        norm_snp = (snp - mean_g) / sd_g
        norm_snp.shape = (num_indivs, 1)

        # Project on the PCs
        pcs += sp.dot(norm_snp, pc_weights)
        num_snps_used += 1

    return {'num_snps_used': num_snps_used, 'num_nt_issues': num_nt_issues, 'pcs': pcs}

def test_assignment4():
    ##Example without noise
    x_train = sp.array([[ 0,  0,  1 , 1],[ 0,  1,  0, 1]])
    y_train = sp.array([[0, 1, 1, 2]])
    w_est = train_ols(x_train, y_train) 
    w_est_ridge = train_ols(x_train, y_train, llambda = 1)
    assert(sp.all(w_est.T == [[1, 1]])) 
    assert(sp.all(w_est_ridge.T == [[.75, .75]]))
    y_est = apply_ols(w_est,x_train)
    assert(sp.all(y_train == y_est)) 
    print 'No-noise-case tests passed'
	
	##Example with noise
	#Data generation
    w_true = 4
    X_train = sp.arange(10)
    X_train = X_train[None,:]
    Y_train = w_true * X_train + sp.random.normal(0,2,X_train.shape)
    #Regression 
    w_est = train_ols(X_train, Y_train) 
    Y_est = apply_ols(w_est,X_train)
    #Plot result
    pl.figure()
    pl.plot(X_train.T, Y_train.T, '+', label = 'Train Data')
    pl.plot(X_train.T, Y_est.T, label = 'Estimated regression')
    pl.xlabel('x')
    pl.ylabel('y')
    pl.legend(loc = 'lower right')

def est_condprob2(data, val, given):
    """Calculate the probability of P(X|Y,Z)

    est_condprob2(data, 'A', ['M', 'LC'])"""

    if not isinstance(given, list):
        raise IndexError("Given must be a list or tuple of givens")
    elif len(given) != 2:
        raise IndexError("I need multiple givens!  Give me more...give me more!")

    gcols = []
    for g in given:
        if g in ['M', 'F']:
            gcols.append(1)
        elif g in ['LC', 'SC', 'T']:
            gcols.append(2)
        elif g in ['A', 'B', 'C']:
            gcols.append(0)

    if val in ['M', 'F']:
        vcol = 1
    elif val in ['LC', 'SC', 'T']:
        vcol = 2
    elif val in ['A', 'B', 'C']:
        vcol = 0

    datsize = data.shape[0]
    needed = [val, given[0], given[1]]
    t = sp.where([sp.all(data[i]==needed) for i in range(datsize)])[0]

    t2 = sp.where([sp.all(data[i,1:]==given) for i in range(datsize)])[0]
    
    return float(t.size)/t2.size

 def testMomentOfInertiaRotatedEllipsoid(self):
     img = mango.zeros(shape=self.imgShape*2, mtype="tomo", origin=(0,0,0))
     img.md.setVoxelSize((1,1,1))
     img.md.setVoxelSizeUnit("mm")
     c = (sp.array(img.origin) + img.origin + img.shape-1)*0.5
     r = sp.array(img.shape-1)*0.25
     
     mango.data.fill_ellipsoid(img, centre=c, radius=r, fill=512)
     rMatrix = \
         (
             mango.image.rotation_matrix(-25, 2).dot(
             mango.image.rotation_matrix( 10, 1).dot(
             mango.image.rotation_matrix( 45, 0)
             ))
         )
     img = mango.image.affine_transform(img, rMatrix, offset=c-img.origin, interptype=mango.image.InterpolationType.CATMULL_ROM_CUBIC_SPLINE)
     #mango.io.writeDds("tomoMoiRotEllipsoid.nc", img)
     
     pmoi, pmoi_axes, com = mango.image.moment_of_inertia(img)
     rootLogger.info("rmtx = \n%s" % (rMatrix,))
     rootLogger.info("pmoi = \n%s" % (pmoi,))
     rootLogger.info("pmoi_axes = \n%s" % (pmoi_axes,))
     rootLogger.info("c = %s, com = %s" % (c, com))
     self.assertTrue(sp.all(sp.absolute(c - com) <= 1.0e-10))
     self.assertLess(pmoi[0], pmoi[1])
     self.assertLess(pmoi[1], pmoi[2])
     self.assertTrue(sp.all(sp.absolute(pmoi_axes[:,0]-rMatrix[:,2]) <= 1.0e-3))
     self.assertTrue(sp.all(sp.absolute(pmoi_axes[:,1]-rMatrix[:,1]) <= 1.0e-3))
     self.assertTrue(sp.all(sp.absolute(pmoi_axes[:,2]-rMatrix[:,0]) <= 1.0e-3))

 def test_shape_3D(self):
     net = op.network.Cubic(shape=[5, 5, 5])
     assert sp.all(net.shape == [5, 5, 5])
     net = op.network.Cubic(shape=[3, 4, 5])
     assert sp.all(net.shape == [3, 4, 5])
     net = op.network.Cubic(shape=[1, 5, 1])
     assert sp.all(net.shape == [1, 5, 1])

    def check( x ):
        y = sl.canonicalise( x )
        yr = y[0,:]
        yc = y[:,0]

        assert all( yr == sc.sort( yr ) )
        assert all( yc == sc.sort( yc ) )

 def test_map_pores_reverse(self):
     Ps = self.net2.Ps[:5]
     b = self.geo21.map_pores(pores=Ps, origin=self.net2)
     assert sp.all(b == [0, 1, 2, 3, 4])
     Ps = self.net2.Ps[-5:]
     b = self.geo22.map_pores(pores=Ps, origin=self.net2)
     assert sp.all(b == [4, 5, 6, 7, 8])

def effectiveHeight(surf1, surf2, plasma, t, lim1=.88, lim2=.92):
    """ calculate the effective height of a view through the scrape-off-layer"""

    segments = viewPoints(surf1,surf2,plasma,t,lim1=lim1,lim2=lim2,fillorder=False)
    output = scipy.zeros((len(segments),))
    for i in xrange(len(segments)):
        area = 0
        if not scipy.all(segments[i] == []):
            
            inlen = geometry.pts2Vec(segments[i][0][0][0],segments[i][0][0][1])
            outlen = geometry.pts2Vec(segments[i][1][0][0],segments[i][1][0][1])

 
            for j in xrange(len(segments[i])): # loop over in vs out

                temp = []
                for k in segments[i][j]: # loop over number of 'rays'
                    for l in k:
                        temp += [l.x()[[0,2]]]
                temp = scipy.array(temp)
        #delaunay
                for k in xrange((len(temp)-2)/2):
                    y = temp[[0,1,2*(k+1),2*(k+1)+1]]
                    if not(scipy.all([y[0] == y[1],y[2] == y[3]])):
                        tri = scipy.spatial.Delaunay(y) #divide into a lower and upper section
               
                #plt.triplot(y[:,0],y[:,1],tri.vertices.copy())
                        for l in tri.points[tri.vertices]:
                            area += calcArea(l)


                output[i] = area/(inlen.s + outlen.s)               
    #plt.show()
    return output

 def test_late_pore_filling(self):
     self.phase['pore.pc_star'] = 1000
     self.phys.add_model(propname='pore.nwp_saturation',
                         model=pm.multiphase.late_filling,
                         Pc_star='pore.pc_star',
                         pressure='pore.pressure')
     assert sp.all(self.phase['pore.nwp_saturation'] < 1.0)
     assert sp.all(self.phase['pore.nwp_saturation'] > 0.0)

def isequal(A, B):
    
    if A is None or B is None:
        return False

    if not sp.all(A.shape == B.shape):
        return False
    else:
        return sp.all(A == B)

 def isDuplicateImage(self, md, mdOther):
     return \
         (
             sp.all((md.voxelSize - mdOther.voxelSize) < 1.0e-5)
             and
             sp.all(md.voxelUnit == mdOther.voxelUnit)
             and
             sp.all(md.volumeSize == mdOther.volumeSize)
         )

 def test_channels(self):
     self.Data.data[...] = 10
     self.DataSub.data[...] = 10 + rand.rand(*self.Data.data.shape)
     self.DataSub.data[:,:,:,123] = 10 + 6*rand.rand(
             *self.Data.data.shape[:-1])
     self.DataSub.data[:,:,:,54] += 10
     reflag.flag(self.Data, self.DataSub, thres=2.0, max_noise_factor=3)
     self.assertTrue(sp.all(self.Data.data.mask[:,:,:,123]))
     self.assertFalse(sp.all(self.Data.data.mask[:,:,:,54]))

    def test_times_limited_by_t_start_and_t_stop(self):
        t_start = 10 * pq.s
        t_stop = 20 * pq.s

        st = self.invoke_gen_func(
            self.defaultRate, t_start=t_start, t_stop=t_stop)
        self.assertTrue(sp.all(t_start < st))
        self.assertTrue(sp.all(st <= t_stop))
        self.assertEqual(t_start, st.t_start)
        self.assertEqual(t_stop, st.t_stop)

def _compare_hdf5(expected, actual):
    for k in expected:
        if isinstance(expected[k], h5py._hl.group.Group):
            for l in expected[k]:
                assert sp.all(expected[k][l][:] == actual[k][l][:])
        else:
            if k in ['psi']:
                assert sp.all(expected[k][:].astype('str') == actual[k][:].astype('str'))
            else:
                assert sp.all(expected[k][:] == actual[k][:])