void GEPPGrowth( Matrix<T>& A, Int n )
{
    DEBUG_ONLY(CSE cse("GEPPGrowth"))
    Identity( A, n, n );
    if( n <= 1 )
        return;

    // Set the last column to all ones
    auto aLast = A( IR(0,n), IR(n-1,n) );
    Fill( aLast, T(1) );

    // Set the subdiagonals to -1
    for( Int j=1; j<n; ++j )
        FillDiagonal( A, T(-1), -j );
}

    /**
     * @brief SetScaleMatrix set the matrix as a scaling matrix.
     * @param x_scale
     * @param y_scale
     */
    void SetScaleMatrix(float x_scale, float y_scale)
    {
        if(x_scale <= 0.0f) {
            x_scale = 1.0f;
        }

        if(y_scale <= 0.0f) {
            y_scale = 1.0f;
        }

        Identity();

        data[0] = x_scale;
        data[4] = y_scale;
    }

inline void
ExplicitLQHelper( Matrix<Complex<Real> >& L, Matrix<Complex<Real> >& A )
{
    Matrix<Complex<Real> > t;
    LQ( A, t );
    L = A;
    MakeTrapezoidal( LEFT, LOWER, 0, L );

    // TODO: Replace this with an in-place expansion of Q
    Matrix<Complex<Real> > Q;
    Identity( A.Height(), A.Width(), Q );
    ApplyPackedReflectors
    ( RIGHT, UPPER, HORIZONTAL, BACKWARD, UNCONJUGATED, 0, A, t, Q );
    A = Q;
}

inline void
NewtonSchulzStep( const Matrix<F>& X, Matrix<F>& XTmp, Matrix<F>& XNew )
{
#ifndef RELEASE
    CallStackEntry entry("sign::NewtonSchulzStep");
#endif
    typedef BASE(F) Real;
    const Int n = X.Height();
 
    // XTmp := 3I - X^2
    Identity( XTmp, n, n );
    Gemm( NORMAL, NORMAL, Real(-1), X, X, Real(3), XTmp );

    // XNew := 1/2 X XTmp
    Gemm( NORMAL, NORMAL, Real(1)/Real(2), X, XTmp, XNew );
}

Expression *IdentityExp::optimize(int result)
{
    //printf("IdentityExp::optimize(result = %d) %s\n", result, toChars());
    e1 = e1->optimize(WANTvalue | (result & WANTinterpret));
    e2 = e2->optimize(WANTvalue | (result & WANTinterpret));
    Expression *e = this;

    if ((this->e1->isConst()     && this->e2->isConst()) ||
        (this->e1->op == TOKnull && this->e2->op == TOKnull))
    {
        e = Identity(op, type, this->e1, this->e2);
        if (e == EXP_CANT_INTERPRET)
            e = this;
    }
    return e;
}

    static Matrix4<T> Rotate(T radians, const vec3& axis) {
        T s = std::sin(radians);
        T c = std::cos(radians);

        Matrix4 m = Identity();
        m.x.x = c + (1 - c) * axis.x * axis.x;
        m.x.y = (1 - c) * axis.x * axis.y - axis.z * s;
        m.x.z = (1 - c) * axis.x * axis.z + axis.y * s;
        m.y.x = (1 - c) * axis.x * axis.y + axis.z * s;
        m.y.y = c + (1 - c) * axis.y * axis.y;
        m.y.z = (1 - c) * axis.y * axis.z - axis.x * s;
        m.z.x = (1 - c) * axis.x * axis.z - axis.y * s;
        m.z.y = (1 - c) * axis.y * axis.z + axis.x * s;
        m.z.z = c + (1 - c) * axis.z * axis.z;
        return m;
    }

// Assumes that there is either no exogenous variables or there is a single
// exogenous variable that is constant and equal to one.  The unconditional
// mean is obtained from the reduced form companion matrix.
TDenseMatrix UnconditionalVariance(const TDenseMatrix &A0, const TDenseMatrix &Aplus, bool IsConstant)
{
  int n_lags=NumberLags(A0,Aplus,IsConstant), n_vars=A0.cols;
  TDenseMatrix B=ReducedForm(A0,Aplus);
  if (n_lags == 0) return ConditionalVariance(A0);
  TDenseMatrix C=CompanionMatrix(B,n_lags), V=BlockDiagonalMatrix(A0,n_lags),
    X=V*(Identity(n_vars*n_lags) - C);
  try
    {
      return SubMatrix(Inverse(TransposeMultiply(X,X)),0,n_vars-1,0,n_vars-1);
    }
  catch (dw_exception &e)
    {
      throw dw_exception("UnconditionalMean(): Unconditional mean does not exist");
    }
}

LRESULT ChatCtrl::onReport(WORD /*wNotifyCode*/, WORD /*wID*/, HWND /*hWndCtl*/, BOOL& /*bHandled*/) {
	const OnlineUserPtr ou = client->findUser(Text::fromT(selectedUser));
		
	if(ou) {
			CHARFORMAT2 cf;
			memzero(&cf, sizeof(CHARFORMAT2));
			cf.cbSize = sizeof(cf);
			cf.dwMask = CFM_BACKCOLOR | CFM_COLOR | CFM_BOLD;
			cf.crBackColor = SETTING(BACKGROUND_COLOR);
			cf.crTextColor = SETTING(ERROR_COLOR);

			AppendText(Identity(NULL, 0), Text::toT(client->getCurrentNick()), Text::toT("[" + Util::getShortTimeString() + "] "), Text::toT(WinUtil::getReport(ou->getIdentity(), m_hWnd)) + _T('\n'), cf, false);
		}

	return 0;
}

const ECP::Point& ECP::Add(const Point &P, const Point &Q) const
{
	if (P.identity) return Q;
	if (Q.identity) return P;
	if (GetField().Equal(P.x, Q.x))
		return GetField().Equal(P.y, Q.y) ? Double(P) : Identity();

	FieldElement t = GetField().Subtract(Q.y, P.y);
	t = GetField().Divide(t, GetField().Subtract(Q.x, P.x));
	FieldElement x = GetField().Subtract(GetField().Subtract(GetField().Square(t), P.x), Q.x);
	m_R.y = GetField().Subtract(GetField().Multiply(t, GetField().Subtract(P.x, x)), P.y);

	m_R.x.swap(x);
	m_R.identity = false;
	return m_R;
}

	Matrix4x4 Matrix4x4::Perspective(float fov, float aspect, float zNear, float zFar)
	{
		Matrix4x4 m = Identity();

		float yScale = 1/tan(Mathf::Deg2Rad * fov / 2);
		float xScale = yScale/aspect;

		m.m00 = xScale;
		m.m11 = yScale;
		m.m22 = (zNear + zFar) / (zNear - zFar);
		m.m23 = 2 * zNear * zFar / (zNear - zFar);
		m.m32 = -1.0f;
		m.m33 = 0;

		return m;
	}

const EC2N::Point& EC2N::Double(const Point &P) const
{
	if (P.identity) return P;
	if (!m_field->IsUnit(P.x)) return Identity();

	FieldElement t = m_field->Divide(P.y, P.x);
	m_field->Accumulate(t, P.x);
	m_R.y = m_field->Square(P.x);
	m_R.x = m_field->Square(t);
	m_field->Accumulate(m_R.x, t);
	m_field->Accumulate(m_R.x, m_a);
	m_field->Accumulate(m_R.y, m_field->Multiply(t, m_R.x));
	m_field->Accumulate(m_R.y, m_R.x);

	m_R.identity = false;
	return m_R;
}

inline void
Explicit( Matrix<F>& L, Matrix<F>& A )
{
#ifndef RELEASE
    CallStackEntry cse("lq::Explicit");
#endif
    Matrix<F> t;
    LQ( A, t );
    L = A;
    MakeTriangular( LOWER, L );

    // TODO: Replace this with an in-place expansion of Q
    Matrix<F> Q;
    Identity( Q, A.Height(), A.Width() );
    lq::ApplyQ( RIGHT, NORMAL, A, t, Q );
    A = Q;
}

    void Midpoint_SingleStep(Vector& x,
                                double t,
                                double& tau,
                                Vector& a,
                                double a0,
                                const Matrix& nu,
                                PropensityFunc propFunc,
                                PropensityJacobianFunc propJacFunc,
                                double abs_tol,
                                double rel_tol,
                                int& RXN,
                                Vector& p)
    {
      p = a; 
      // static Vector x1 = x; 
      static Vector a1 = a; 
      // static Vector delx = x; 
/*
      p = PoissonRandom(0.5* a * tau); // better midpoint method
      p = 0.5* a * tau; // original midpoint method
      x1 = x + nu*p; 
*/
      const Vector A = nu * a;
      const Vector E = 0.5* tau * A;
      Vector delx(x.Size(), 0);
      bool converged = false;
      int iterations = 0;
      const Matrix I = Identity(x.Size());
      Vector xPlusDelx(x.Size());
      static Matrix AA = I;
      static Vector BB = E;
      static Vector deldelx = x;
      while (!converged && iterations < 10) {
        xPlusDelx = x+ delx;
        AA = I - 0.5*tau * nu * propJacFunc(xPlusDelx);
        BB = E + 0.5*tau * nu * propFunc(xPlusDelx) - delx;
        SolveGE(AA, deldelx, BB);

        converged = (Norm(deldelx, 2) <= rel_tol * Norm(delx, 2) + abs_tol);
        delx += deldelx;
        ++iterations;
      }
      a1 = propFunc(x + delx); 
      a1 = 0.5*(a + a1)*tau; 
      p = PoissonRandom(a1);
    }

void renderScene(void) {
  int sx = glutGet(GLUT_WINDOW_WIDTH);
  int sy = glutGet(GLUT_WINDOW_HEIGHT);
  //SAspect = Scaling(1, (float)sx / sy, 1);
  VECTOR4D worldRayOrigin, worldRayDir;
  VECTOR4D modelRayOrigin, modelRayDir;
  MATRIX4D InvW;
  multimap<float, CMesh::INTERSECTIONINFO> faces;
  bool fill = false;

  SAspect = Scaling((float)sy / sx, 1, 1);
  W = Identity();

  P = PerspectiveWidthHeightRH(0.5f, 0.5f, 1.0f, 10.0f);
  EC = SAspect * T * Translation(0.0f, 0.0f, -1.0f);

  glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
  BuildRayFromPerspective(EC, mx, my, worldRayOrigin, worldRayDir);
  Inverse(W, InvW);
  modelRayOrigin = InvW * worldRayOrigin;
  modelRayDir = InvW * worldRayDir;
  fill = g_EggCarton.RayCast(modelRayOrigin, modelRayDir, faces);

  glPolygonMode(GL_FRONT_AND_BACK, GL_LINE);
  glBegin(GL_TRIANGLES);
  g_EggCarton.Draw(EC);
  glEnd();
  glPolygonMode(GL_FRONT_AND_BACK, GL_FILL);
  glBegin(GL_TRIANGLES);
  if (fill)
    g_EggCarton.Draw(EC, faces.begin()->second.Face, 1);
  glEnd();

  if (bWireframe)
    glPolygonMode(GL_FRONT_AND_BACK, GL_LINE);
  else
    glPolygonMode(GL_FRONT_AND_BACK, GL_FILL);
  glBegin(GL_TRIANGLES);
  g_Plate.Draw(SAspect * T * Translation(0.0f, 0.0f, -0.1f));
  g_Sphere.Draw(SAspect * T * Translation(0.4f, 0.4f, 0.32f));
  g_Bananas.Draw(SAspect * T * Translation(0.45f, 0.45f, 0.48f));
  g_Flower.Draw(SAspect * T * Translation(-0.5f, 0.5f, 0.8f));
  glEnd();

  glutSwapBuffers();
}

Matrix4 Matrix4::Face(const Vec3f &V0, const Vec3f &V1)
{
    //
    // Rotate about the cross product of the two vectors by the angle between the two vectors
    //
    Vec3f Axis = Vec3f::Cross(V0, V1);
    float Angle = Vec3f::AngleBetween(V0, V1);

    if(Angle == 0.0f || Axis.Length() < 0.0f)
    {
        return Identity();
    }
    else
    {
        return Rotation(Axis, Angle);
    }
}

Matrix3& Matrix3::Orthographic( const Real l, const Real r, const Real b, const Real t, const Real zn, const Real zf )
{
#ifdef __USE_D3DX__
	D3DXMatrixOrthoOffCenterLH( (D3DXMATRIX*)this, l, r, b, t, zn, zf );
#else
	Identity();
		
	Real d = zf - zn;

	e[0] = 2.0f / ( r - l );
	e[5] = 2.0f / ( t - b );
	e[10] = 1.0f / d;
	e[14] = -zn / d;
		
#endif
	return *this;
}

void
NewtonSchulzStep
( const Matrix<Field>& X,
        Matrix<Field>& XTmp,
        Matrix<Field>& XNew )
{
    EL_DEBUG_CSE
    typedef Base<Field> Real;
    const Int n = X.Height();

    // XTmp := 3I - X^2
    Identity( XTmp, n, n );
    Gemm( NORMAL, NORMAL, Real(-1), X, X, Real(3), XTmp );

    // XNew := 1/2 X XTmp
    Gemm( NORMAL, NORMAL, Real(1)/Real(2), X, XTmp, XNew );
}

inline void
ExplicitLQHelper
( DistMatrix<Complex<Real> >& L, DistMatrix<Complex<Real> >& A )
{
    const Grid& g = A.Grid();
    DistMatrix<Complex<Real>,MD,STAR> t( g );
    LQ( A, t );
    L = A;
    MakeTrapezoidal( LEFT, LOWER, 0, L );

    // TODO: Replace this with an in-place expansion of Q
    DistMatrix<Complex<Real> > Q( g );
    Identity( A.Height(), A.Width(), Q );
    ApplyPackedReflectors
    ( RIGHT, UPPER, HORIZONTAL, BACKWARD, UNCONJUGATED, 0, A, t, Q );
    A = Q;
}

// Takes care that there may be multiple values of attribute in valueMap
Identity RecognitionDatabase::Private::findByAttributes(const QString& attribute, const QMap<QString, QString>& valueMap) const
{
    QMap<QString, QString>::const_iterator it = valueMap.find(attribute);

    for (; it != valueMap.end() && it.key() == attribute; ++it)
    {
        foreach (const Identity& identity, identityCache)
        {
            if (identityContains(identity, attribute, it.value()))
            {
                return identity;
            }
        }
    }

    return Identity();
}

  void handleMe(boost::system::error_code err, const Http::Message& response)
  {
#ifndef WT_TARGET_JAVA
    WApplication::instance()->resumeRendering();
#endif

    if (!err && response.status() == 200) {
#ifndef WT_TARGET_JAVA
      Json::ParseError e;
      Json::Object me;
      bool ok = Json::parse(response.body(), me, e);
#else
      Json::Object me;
      try {
	me = (Json::Object)Json::Parser().parse(response.body());
      } catch (Json::ParseError pe) {
      }
      bool ok = me.isNull();
#endif

      if (!ok) {
	LOG_ERROR("could not parse Json: '" << response.body() << "'");
	setError(ERROR_MSG("badjson"));
	authenticated().emit(Identity::Invalid);
      } else {
	std::string id = me.get("id");
	WT_USTRING userName = me.get("name");
	std::string email = me.get("email").orIfNull("");
        bool emailVerified = !me.get("email").isNull();

	authenticated().emit(Identity(service().name(), id, userName,
				      email, emailVerified));
      }
    } else {
      if (!err) {
	LOG_ERROR("user info request returned: " << response.status());
	LOG_ERROR("with: " << response.body());
      } else
	LOG_ERROR("handleMe(): " << err.message());

      setError(ERROR_MSG("badresponse"));

      authenticated().emit(Identity::Invalid);
    }
  }