def pansharpen(imname, inDir, outDir):
    """Create pansharpened images from EO-1 scene directory."""
    pan = gdal.Open(inDir + '/' + imname + '_B' + digits % 1 + '_L1T.TIF')
    bigPan = pan.ReadAsArray()[:-1, :-1]

    bandNums = []
    srcs = []
    dests = []
    driver = gdal.GetDriverByName("GTiff")
    driver.Register()
    for bandNum in range(2, 11):
        tif = gdal.Open(inDir+'/'+imname + '_B' + digits % bandNum + '_L1T.TIF')
        if tif is None:
            print "WARNING: Band %d not found." % bands
            continue
        bandNums.append(bandNum)
        srcs.append(tif)
        dests.append(driver.CreateCopy(outDir + '/' + imname + '_B' +
                                       digits % bandNum + '_L1T.TIF', pan, 0))

    imgs = [np.float64(src.ReadAsArray()[:-1, :-1]) for src in srcs]
    smoothPan = np.kron(sum(imgs) / len(imgs), np.ones([3, 3])) + 1e-9

    for img, dest in zip(imgs, dests):
        newimg = bigPan / smoothPan * np.kron(img, np.ones([3, 3]))
        #newimg[newimg > 255] = 255
        band = dest.GetRasterBand(1)
        band.WriteArray(newimg)
    dest = None
    dests = None

    def act(self, qubits, index=1):
        """Performs the action of a quantum gate on a qubit system.

        Operates in-place on the system "qubits", so the original system is
        changed by interaction with the gate. This avoids violations of the no-
        cloning theorem.

        Args:
            qubits: A system of qubits for the gate to act on.
            index: Starting index of the first qubit in the system for the gate
            to act one. E.g. if the gate acts on two qubits and index = 3, then
            the gate would act on the 3rd and 4th qubits in the system (where
            qubits are indexed starting at one).

        Raises:
            ValueError if there is a dimension mismatch between the gate and the
            qubits to be operated on.
            RuntimeError if the matrix is not unitary.
        """
        if index <= 0 or index + self.n() - 1 > qubits.n():
            raise ValueError('Dimension mismatch with gate and qubit system')
        if not is_unitary(self.__matrix):
            raise RuntimeError('Non-unitary matrix')

        # construct a matrix to operate only on the desired qubits
        G = copy(self.__matrix)
        if index > 1:
            G = numpy.kron(eye(2**(index - 1)), G)
        if index + self.n() - 1 < qubits.n():
            G = numpy.kron(G, eye(2**(qubits.n() - (index + self.n() -1))))

        qubits._QubitSystem__coeffs = numpy.dot(G, qubits._QubitSystem__coeffs)

	def updateParameters(self, articlePicked, click,  userID):	
		self.counter +=1
		self.Wlong = vectorize(self.W)
		featureDimension = len(articlePicked.featureVector)
		T_X = vectorize(np.outer(articlePicked.featureVector, self.W.T[userID])) 
		self.A += np.outer(T_X, T_X)	
		self.b += click*T_X
		self.AInv = np.linalg.inv(self.A)
		self.UserTheta = matrixize(np.dot(self.AInv, self.b), len(articlePicked.featureVector)) 

		Xi_Matirx = np.zeros(shape = (featureDimension, self.userNum))
		Xi_Matirx.T[userID] = articlePicked.featureVector
		W_X = vectorize( np.dot(np.transpose(self.UserTheta), Xi_Matirx))
		self.batchGradient +=evaluateGradient(W_X, click, self.Wlong, self.lambda_, self.regu  )

		if self.counter%self.windowSize ==0:
			self.Wlong -= 1/(float(self.counter/self.windowSize)+1)*self.batchGradient
			self.W = matrixize(self.Wlong, self.userNum)
			self.W = normalize(self.W, axis=0, norm='l1')
			#print 'SVD', self.W
			self.batchGradient = np.zeros(self.userNum*self.userNum)
			# Use Ridge regression to fit W
		'''
		plt.pcolor(self.W_b)
		plt.colorbar
		plt.show()
		'''
		if self.W.T[userID].any() <0 or self.W.T[userID].any()>1:
			print self.W.T[userID]

		self.CoTheta = np.dot(self.UserTheta, self.W)
		self.BigW = np.kron(np.transpose(self.W), np.identity(n=len(articlePicked.featureVector)))
		self.CCA = np.dot(np.dot(self.BigW , self.AInv), np.transpose(self.BigW))
		self.BigTheta = np.kron(np.identity(n=self.userNum) , self.UserTheta)

def expand(A):
    M,N = A.shape
    t = np.kron(A.flatten(),np.ones(N))
    u = np.triu(np.ones((N,N))).flatten()
    v = np.kron(np.ones(M),u)
    w  = t *  v
    return(w.reshape(M,N,N).swapaxes(1,2))

    def test_vcomp_3(self):
        # Test a model with vcomp but no other random effects, using formulas.

        import scipy

        v = scipy.__version__.split(".")[1]
        v = int(v)
        if v < 16:
            return

        np.random.seed(4279)
        x1 = np.random.normal(size=400)
        groups = np.kron(np.arange(100), np.ones(4))
        slopes = np.random.normal(size=100)
        slopes = np.kron(slopes, np.ones(4)) * x1
        y = slopes + np.random.normal(size=400)
        vc_fml = {"a": "0 + x1"}
        df = pd.DataFrame({"y": y, "x1": x1, "groups": groups})

        model = MixedLM.from_formula("y ~ 1", groups="groups", vc_formula=vc_fml, data=df)
        result = model.fit()
        result.summary()

        assert_allclose(result.resid.iloc[0:4], np.r_[-1.180753, 0.279966, 0.578576, -0.667916], rtol=1e-3)
        assert_allclose(result.fittedvalues.iloc[0:4], np.r_[-0.101549, 0.028613, -0.224621, -0.126295], rtol=1e-3)

def test_symm_algorithm_equivalence():
    """Test different stabilization methods in the computation of modes,
    in the presence and/or absence of the discrete symmetries."""
    np.random.seed(400)
    for n in (12, 20, 40):
        for sym in kwant.rmt.sym_list:
            # Random onsite and hopping matrices in symmetry class
            h_cell = kwant.rmt.gaussian(n, sym)
            # Hopping is an offdiagonal block of a Hamiltonian. We rescale it
            # to ensure that there are modes at the Fermi level.
            h_hop = 10 * kwant.rmt.gaussian(2*n, sym)[:n, n:]

            if kwant.rmt.p(sym):
                p_mat = np.array(kwant.rmt.h_p_matrix[sym])
                p_mat = np.kron(np.identity(n // len(p_mat)), p_mat)
            else:
                p_mat = None

            if kwant.rmt.t(sym):
                t_mat = np.array(kwant.rmt.h_t_matrix[sym])
                t_mat = np.kron(np.identity(n // len(t_mat)), t_mat)
            else:
                t_mat = None

            if kwant.rmt.c(sym):
                c_mat = np.kron(np.identity(n // 2), np.diag([1, -1]))
            else:
                c_mat = None

            check_equivalence(h_cell, h_hop, n, sym=sym, particle_hole=p_mat,
                              chiral=c_mat, time_reversal=t_mat)

    def _make_pairs(self, i, j):
        """
        Create arrays containing all unique ordered pairs of i, j.

        The arrays i and j must be one-dimensional containing non-negative
        integers.
        """

        mat = np.zeros((len(i) * len(j), 2), dtype=np.int32)

        # Create the pairs and order them
        f = np.ones(len(j))
        mat[:, 0] = np.kron(f, i).astype(np.int32)
        f = np.ones(len(i))
        mat[:, 1] = np.kron(j, f).astype(np.int32)
        mat.sort(1)

        # Remove repeated rows
        try:
            dtype = np.dtype((np.void, mat.dtype.itemsize * mat.shape[1]))
            bmat = np.ascontiguousarray(mat).view(dtype)
            _, idx = np.unique(bmat, return_index=True)
        except TypeError:
            # workaround for old numpy that can't call unique with complex
            # dtypes
            rs = np.random.RandomState(4234)
            bmat = np.dot(mat, rs.uniform(size=mat.shape[1]))
            _, idx = np.unique(bmat, return_index=True)
        mat = mat[idx, :]

        return mat[:, 0], mat[:, 1]

def getU(J,hx,hz,t):
    """
    Time evolution operator acting on 2 spins.

    Parameters
    ----------
    J : float
        Pair-wise interaction energy.
    hx : float
        Magnetic energy of each spin with dipole moment mu in field B.
    t : float
        Timestep of each iteration.

    Returns
    -------
    U : (2,2,2,2) ndarray
        Non-unitary time evolution operator.        
    """
    I = s(0)
    X = s(1)
    Z = s(3)
    hamiltonian = -J*np.kron(Z,Z)\
           -(np.kron(X,I)+np.kron(I,X))*hx*0.5\
           -(np.kron(Z,I)+np.kron(I,Z))*hz*0.5
    U = expm(-hamiltonian*t)
    return np.reshape(U,(2,2,2,2))

def pauli(paulis,positions,N):
    """
    N-qubit Pauli operator given paulis and positions.

    Parameters
    ----------
    paulis : list
        List of integers denoting type of Pauli matrix at each site.
    positions : list
        List of positions for each corresponding element in paulis.
    N : int
        Total number of sites.

    Returns
    -------
    pauli : (2**N,2**N) ndarray
        Pauli operator as 2^N by 2^N matrix.
    """
    mat = 1+0j
    identity = s(0)
    for i in xrange(N):
        if i in positions:
            mat = np.kron(mat,s(paulis[positions.index(i)]))
        else:
            mat = np.kron(mat,identity)
    return mat

def test_qubit_order_to_wavefunction_order_matches_np_kron(scheduler):
    simulator = cg.XmonSimulator()
    zero = [1, 0]
    one = [0, 1]

    result = simulate(simulator,
                      cirq.Circuit.from_ops(cirq.X(Q1)),
                      scheduler,
                      qubit_order=[Q1, Q2])
    assert cirq.allclose_up_to_global_phase(
        result.final_state, np.kron(one, zero))

    result = simulate(simulator,
                      cirq.Circuit.from_ops(cirq.X(Q1)),
                      scheduler,
                      qubit_order=[Q2, Q1])
    assert cirq.allclose_up_to_global_phase(
        result.final_state, np.kron(zero, one))

    result = simulate(simulator,
                      cirq.Circuit.from_ops(cirq.X(Q1)),
                      scheduler,
                      qubit_order=cirq.QubitOrder.sorted_by(repr))
    assert cirq.allclose_up_to_global_phase(
        result.final_state, np.array(one))

    result = simulate(simulator,
                      cirq.Circuit.from_ops(cirq.X(Q1), cirq.Z(Q2)),
                      scheduler,
                      qubit_order=cirq.QubitOrder.sorted_by(repr))
    assert cirq.allclose_up_to_global_phase(
        result.final_state, np.kron(one, zero))

    def test_summary(self):
        # smoke test
        np.random.seed(34234)
        time = 50 * np.random.uniform(size=200)
        status = np.random.randint(0, 2, 200).astype(np.float64)
        exog = np.random.normal(size=(200,4))

        mod = PHReg(time, exog, status)
        rslt = mod.fit()
        smry = rslt.summary()

        strata = np.kron(np.arange(50), np.ones(4))
        mod = PHReg(time, exog, status, strata=strata)
        rslt = mod.fit()
        smry = rslt.summary()
        msg = "3 strata dropped for having no events"
        assert_(msg in str(smry))

        groups = np.kron(np.arange(25), np.ones(8))
        mod = PHReg(time, exog, status)
        rslt = mod.fit(groups=groups)
        smry = rslt.summary()

        entry = np.random.uniform(0.1, 0.8, 200) * time
        mod = PHReg(time, exog, status, entry=entry)
        rslt = mod.fit()
        smry = rslt.summary()
        msg = "200 observations have positive entry times"
        assert_(msg in str(smry))

    def liouvillian(self, chi=None):
        sys_hamiltonian = (
            np.kron(self.electronic_hamiltonian(), np.eye(self.vib_basis_size))
            + np.kron(np.eye(3), self.vibrational_hamiltonian())
            + self.el_vib_interaction_hamiltonian()
        )
        L = os.super_operator(sys_hamiltonian, self.lead_operators() + self.vib_damping_operators())

        if chi:
            L_jump = self.jump_liouvillian()
            for i, row in enumerate(L_jump):
                for j, v in enumerate(row):
                    if v != 0:
                        L[i, j] *= np.exp(chi)

        if self.remove_elements:
            L = np.delete(L, self.element_indices_to_remove, 0)
            L = np.delete(L, self.element_indices_to_remove, 1)

            # check physicality of Liouvillian, only need to check that columns related to populations add to 1 (maybe there is a rule for coherences too?)
            dv_pops = np.eye(3 * self.vib_basis_size).flatten()
            dv_pops = np.delete(dv_pops, self.element_indices_to_remove, 0)
            trans_L = L.T
            for i, el in enumerate(dv_pops):
                if el == 1:
                    sum = np.sum(trans_L[i])
                    if not -1.0e-6 < sum < 1.0e-6:
                        print "Liouvillian not physical!"
                        if chi != 0:
                            print "but chi is non zero"
                        print sum

        return L

def gen_crossed_logit_pandas(nc, cs, s1, s2):

    np.random.seed(3799)

    a = np.kron(np.arange(nc), np.ones(cs))
    b = np.kron(np.ones(cs), np.arange(nc))
    fe = np.ones(nc * cs)

    vc = np.zeros(nc * cs)
    for i in np.unique(a):
        ii = np.flatnonzero(a == i)
        vc[ii] += s1*np.random.normal()
    for i in np.unique(b):
        ii = np.flatnonzero(b == i)
        vc[ii] += s2*np.random.normal()

    lp = -0.5 * fe + vc
    pr = 1 / (1 + np.exp(-lp))
    y = 1*(np.random.uniform(size=nc*cs) < pr)

    ident = np.zeros(2*nc, dtype=np.int)
    ident[nc:] = 1

    df = pd.DataFrame({"fe": fe, "a": a, "b": b, "y": y})

    return df

def generate_inequalities_constraints_mat(N, nx, nu, xmin, xmax, umin, umax):
    """
    generate matrices of inequalities constrints

    return G, h
    """
    G = np.zeros((0, (nx + nu) * N))
    h = np.zeros((0, 1))
    if umax is not None:
        tG = np.hstack([np.eye(N * nu), np.zeros((N * nu, nx * N))])
        th = np.kron(np.ones((N * nu, 1)), umax)
        G = np.vstack([G, tG])
        h = np.vstack([h, th])

    if umin is not None:
        tG = np.hstack([np.eye(N * nu) * -1.0, np.zeros((N * nu, nx * N))])
        th = np.kron(np.ones((N, 1)), umin * -1.0)
        G = np.vstack([G, tG])
        h = np.vstack([h, th])

    if xmax is not None:
        tG = np.hstack([np.zeros((N * nx, nu * N)), np.eye(N * nx)])
        th = np.kron(np.ones((N, 1)), xmax)
        G = np.vstack([G, tG])
        h = np.vstack([h, th])

    if xmin is not None:
        tG = np.hstack([np.zeros((N * nx, nu * N)), np.eye(N * nx) * -1.0])
        th = np.kron(np.ones((N, 1)), xmin * -1.0)
        G = np.vstack([G, tG])
        h = np.vstack([h, th])

    return G, h

    def txest_vcomp_1(self):
        """
        Fit the same model using constrained random effects and variance components.
        """

        np.random.seed(4279)
        exog = np.random.normal(size=(400, 1))
        exog_re = np.random.normal(size=(400, 2))
        groups = np.kron(np.arange(100), np.ones(4))
        slopes = np.random.normal(size=(100, 2))
        slopes[:, 1] *= 2
        slopes = np.kron(slopes, np.ones((4, 1))) * exog_re
        errors = slopes.sum(1) + np.random.normal(size=400)
        endog = exog.sum(1) + errors

        free = MixedLMParams(1, 2, 0)
        free.fe_params = np.ones(1)
        free.cov_re = np.eye(2)
        free.vcomp = np.zeros(0)

        model1 = MixedLM(endog, exog, groups, exog_re=exog_re)
        result1 = model1.fit(free=free)

        exog_vc = {"a": {}, "b": {}}
        for k,group in enumerate(model1.group_labels):
            ix = model1.row_indices[group]
            exog_vc["a"][group] = exog_re[ix, 0:1]
            exog_vc["b"][group] = exog_re[ix, 1:2]
        model2 = MixedLM(endog, exog, groups, exog_vc=exog_vc)
        result2 = model2.fit()
        result2.summary()

        assert_allclose(result1.fe_params, result2.fe_params, atol=1e-4)
        assert_allclose(np.diag(result1.cov_re), result2.vcomp, atol=1e-2, rtol=1e-4)
        assert_allclose(result1.bse[[0, 1, 3]], result2.bse, atol=1e-2, rtol=1e-2)

def SENSOR(d1=0.1, d2=0.1, d3=0.1):
    L1, L2, L3 = Ls(d1, d2, d3)
    LL1 = np.dot(PDIAG, np.kron(L1, L1))
    LL2 = np.dot(PDIAG, np.kron(L2, L2))
    LL3 = np.dot(PDIAG, np.kron(L3, L3))
    SENSOR = np.r_[LL1[[0, 4, 8], :], LL2[[0, 4, 8], :], LL3[[0, 4, 8], :]]
    return SENSOR

    def test_formulas(self):
        np.random.seed(2410)
        exog = np.random.normal(size=(300, 4))
        exog_re = np.random.normal(size=300)
        groups = np.kron(np.arange(100), [1, 1, 1])
        g_errors = exog_re * np.kron(np.random.normal(size=100),
                                     [1, 1, 1])
        endog = exog.sum(1) + g_errors + np.random.normal(size=300)

        mod1 = MixedLM(endog, exog, groups, exog_re)
        # test the names
        assert_(mod1.data.xnames == ["x1", "x2", "x3", "x4"])
        assert_(mod1.data.exog_re_names == ["x_re1"])
        assert_(mod1.data.exog_re_names_full == ["x_re1 RE"])
        rslt1 = mod1.fit()

        # Fit with a formula, passing groups as the actual values.
        df = pd.DataFrame({"endog": endog})
        for k in range(exog.shape[1]):
            df["exog%d" % k] = exog[:, k]
        df["exog_re"] = exog_re
        fml = "endog ~ 0 + exog0 + exog1 + exog2 + exog3"
        re_fml = "0 + exog_re"
        mod2 = MixedLM.from_formula(fml, df, re_formula=re_fml,
                                    groups=groups)

        assert_(mod2.data.xnames == ["exog0", "exog1", "exog2", "exog3"])
        assert_(mod2.data.exog_re_names == ["exog_re"])
        assert_(mod2.data.exog_re_names_full == ["exog_re RE"])

        rslt2 = mod2.fit()
        assert_almost_equal(rslt1.params, rslt2.params)

        # Fit with a formula, passing groups as the variable name.
        df["groups"] = groups
        mod3 = MixedLM.from_formula(fml, df, re_formula=re_fml,
                                    groups="groups")
        assert_(mod3.data.xnames == ["exog0", "exog1", "exog2", "exog3"])
        assert_(mod3.data.exog_re_names == ["exog_re"])
        assert_(mod3.data.exog_re_names_full == ["exog_re RE"])

        rslt3 = mod3.fit(start_params=rslt2.params)
        assert_allclose(rslt1.params, rslt3.params, rtol=1e-4)

        # Check default variance structure with non-formula model
        # creation, also use different exog_re that produces a zero
        # estimated variance parameter.
        exog_re = np.ones(len(endog), dtype=np.float64)
        mod4 = MixedLM(endog, exog, groups, exog_re)
        with warnings.catch_warnings():
            warnings.simplefilter("ignore")
            rslt4 = mod4.fit()
        from statsmodels.formula.api import mixedlm
        mod5 = mixedlm(fml, df, groups="groups")
        assert_(mod5.data.exog_re_names == ["groups"])
        assert_(mod5.data.exog_re_names_full == ["groups RE"])
        with warnings.catch_warnings():
            warnings.simplefilter("ignore")
            rslt5 = mod5.fit()
        assert_almost_equal(rslt4.params, rslt5.params)

def find_activity_intervals(C, Npeaks=5, tB=-3, tA=10, thres=0.3):
    # todo todocument
    import peakutils
    K, T = np.shape(C)
    L = []
    for i in range(K):
        if np.sum(np.abs(np.diff(C[i, :]))) == 0:
            L.append([])
            print('empyty component at:' + str(i))
            continue
        indexes = peakutils.indexes(C[i, :], thres=thres)
        srt_ind = indexes[np.argsort(C[i, indexes])][::-1]
        srt_ind = srt_ind[:Npeaks]
        L.append(srt_ind)

    LOC = []
    for i in range(K):
        if len(L[i]) > 0:
            interval = np.kron(L[i], np.ones(int(np.round(tA - tB)), dtype=int)) + \
                np.kron(np.ones(len(L[i]), dtype=int), np.arange(tB, tA))
            interval[interval < 0] = 0
            interval[interval > T - 1] = T - 1
            LOC.append(np.array(list(set(interval))))
        else:
            LOC.append(None)

    return LOC

def super_operator(H, jump_operators):
    #print "Constructing super-operator..."
    I = np.eye(H.shape[0], H.shape[1], dtype='complex')
    L = -1.j * (np.kron(H, I) - np.kron(I, H))
    if jump_operators:
        L += incoherent_super_operator(jump_operators)
    return L

def test_measurement_state_update(num_prefix_qubits):
    np.random.seed(3)
    with xmon_stepper.Stepper(num_qubits=3,
                              num_prefix_qubits=num_prefix_qubits,
                              min_qubits_before_shard=0) as s:
        # 1/sqrt(2)(I+iX) gate.
        for i in range(3):
            s.simulate_w(i, -0.5, 0)
        # Check state before measurements.
        single_qubit_state = np.array([1, 1j]) / np.sqrt(2)
        two_qubit_state = np.kron(single_qubit_state, single_qubit_state)
        expected = np.kron(two_qubit_state, single_qubit_state).flatten()
        np.testing.assert_almost_equal(expected, s.current_state)
        assert not s.simulate_measurement(0)
        # Check state after collapse of first qubit state.
        expected = np.kron(two_qubit_state, np.array([1, 0])).flatten()
        np.testing.assert_almost_equal(expected, s.current_state)
        assert not s.simulate_measurement(1)
        # Check state after collapse of second qubit state.
        expected = np.kron(single_qubit_state,
                           np.array([1, 0, 0, 0])).flatten()
        np.testing.assert_almost_equal(expected, s.current_state, decimal=6)
        assert s.simulate_measurement(2)
        expected = np.array([0, 0, 0, 0, 1j, 0, 0, 0])
        # Check final state after collapse of third qubit state.
        np.testing.assert_almost_equal(expected, s.current_state)