C++ InnerProduct函数代码示例

bool CAR_DKW_o::SetParameters_InitializeABOmega()
{
	if (!CAR_DKW::SetParameters_InitializeABOmega())
		return false; 

	TDenseMatrix lambda1 = SIGMA_inverse * SIGMAlambda1;
        TDenseVector KAPPAtheta = KAPPA * theta;
        TDenseVector theta_Q;

	aI_Q.Zeros(dataP->MATgrid_options.Dimension());
        bI_Q.Zeros(dataP->MATgrid_options.Dimension(), Nfac);
        
        for (int i=0; i<dataP->MATgrid_options.Dimension(); i++)
        {
                double MAT = dataP->MATgrid_options(i);
                TDenseVector temp_ay;
                TDenseMatrix temp_by;
                if (!YieldFacLoad(temp_ay, temp_by, KAPPA_rn, Inv_KAPPA_rn, Inv_Kron_KAPPA_rn, SIGMA, KAPPAtheta, rho0, rho1, lambda0,TDenseVector(1,MAT)))
                        return false;
                TDenseVector temp_by_vector = temp_by.RowVector(0); // -MAT * temp_by.RowVector(0);

                theta_Q = Multiply(Inv_KAPPA_rn, KAPPAtheta-SIGMA*lambda0+MultiplyTranspose(SIGMA,SIGMA)*temp_by_vector);

                double rho0_Q = rho0_pi - InnerProduct(lambda0, sigq)+InnerProduct(sigq, TransposeMultiply(SIGMA,temp_by_vector));
                TDenseVector rho1_Q = rho1_pi - TransposeMultiply(lambda1, sigq);

                double temp_aI_Q;
                TDenseVector temp_bI_Q;
                InfExpFacLoad(temp_aI_Q, temp_bI_Q, KAPPA_rn, Inv_KAPPA_rn, Inv_Kron_KAPPA_rn, SIGMA, theta_Q, sigq, sigqx, rho0_Q, rho1_Q, MAT);

                aI_Q(i) = temp_aI_Q;
                bI_Q.InsertRowMatrix(i, 0, temp_bI_Q);
        }
	return true; 
}

void CompMod(zz_pX& x, const zz_pX& g, const zz_pXArgument& A, 
             const zz_pXModulus& F)
{
   if (deg(g) <= 0) {
      x = g;
      return;
   }


   zz_pX s, t;
   vec_zz_p scratch(INIT_SIZE, F.n);

   long m = A.H.length() - 1;
   long l = ((g.rep.length()+m-1)/m) - 1;

   zz_pXMultiplier M;
   build(M, A.H[m], F);

   InnerProduct(t, g.rep, l*m, l*m + m - 1, A.H, F.n, scratch);
   for (long i = l-1; i >= 0; i--) {
      InnerProduct(s, g.rep, i*m, i*m + m - 1, A.H, F.n, scratch);
      MulMod(t, t, M, F);
      add(t, t, s);
   }

   x = t;
}

void PlainUpdateMap(vec_zz_p& xx, const vec_zz_p& a, 
                    const zz_pX& b, const zz_pX& f)
{
   long n = deg(f);
   long i, m;

   if (IsZero(b)) {
      xx.SetLength(0);
      return;
   }

   m = n-1 - deg(b);

   vec_zz_p x(INIT_SIZE, n);

   for (i = 0; i <= m; i++)
      InnerProduct(x[i], a, b.rep, i);

   if (deg(b) != 0) {
      zz_pX c(INIT_SIZE, n);
      LeftShift(c, b, m);

      for (i = m+1; i < n; i++) {
         MulByXMod(c, c, f);
         InnerProduct(x[i], a, c.rep);
      }
   }

   xx = x;
}

void CompMod(ZZ_pX& x, const ZZ_pX& g, const ZZ_pXArgument& A, 
             const ZZ_pXModulus& F)
{
   if (deg(g) <= 0) {
      x = g;
      return;
   }


   ZZ_pX s, t;
   ZZVec scratch(F.n, ZZ_p::ExtendedModulusSize());

   long m = A.H.length() - 1;
   long l = ((g.rep.length()+m-1)/m) - 1;

   ZZ_pXMultiplier M;
   build(M, A.H[m], F);

   InnerProduct(t, g.rep, l*m, l*m + m - 1, A.H, F.n, scratch);
   for (long i = l-1; i >= 0; i--) {
      InnerProduct(s, g.rep, i*m, i*m + m - 1, A.H, F.n, scratch);
      MulMod(t, t, M, F);
      add(t, t, s);
   }

   x = t;
}

double DoubleDeriv1 (edge *TheEdge, position *ThePosition, line *TheLine,
		     int index_x, int index_y)
{
  int index;
  double dv[3], v[3], newv[3], hat[3], temp1;

  if (!TheEdge || !ThePosition || !TheLine) return (0.0);

  if (index_y < index_x) {
    index = index_x; index_x = index_y; index_y = index;
  }

  ComputeV (v, *TheLine);
  RotateQuaternion (newv, ThePosition->q, v);
  temp1 = InnerProduct (newv, TheEdge->m);

  if ((index_x == THETA) && (index_y == THETA) ||
      (index_x == OMEGA) && (index_y == OMEGA))
    return (-2.0 * sqr(temp1));

  if ((index_x == THETA) && (index_y == OMEGA))
    return (0.0);

  if (((index_x == THETA) || (index_x == OMEGA)) &&
      ((index_y == ALPHA) || (index_y == BETA) || (index_y == GAMMA))) {

    if (index_x == THETA) {
      hat[0] = 0.0; hat[1] = -1.0; hat[2] = 0.0;
    } else {
      hat[0] = 1.0; hat[1] =  0.0; hat[2] = 0.0;
    }

    RotateQuaternion (v, TheLine->q, hat);
    RotateQuaternion (newv, ThePosition->q, v);

    DerivR (dv, newv, IndexToAxis (index_y));

    return (2.0 * temp1 * InnerProduct (dv, TheEdge->m));
  }

  if ( ((index_x == ALPHA) || (index_x == BETA) || (index_x == GAMMA)) &&
       ((index_y == ALPHA) || (index_y == BETA) || (index_y == GAMMA)) ) {
    
    ComputeV (v, *TheLine);
    RotateQuaternion (newv, ThePosition->q, v);

    DoubleDerivR (dv, newv, IndexToAxis (index_x), IndexToAxis (index_y));

    return (2.0 * temp1 * InnerProduct (dv, TheEdge->m));
  }

  pexit ("Bad Parameters to DoubleDeriv1");

  return 0.0;
}

static
void ComputeGS(mat_ZZ& B, xdouble **B1, xdouble **mu, xdouble *b, 
               xdouble *c, long k, xdouble bound, long st, xdouble *buf)
{
   long n = B.NumCols();
   long i, j;
   xdouble s, t1, y, t;
   ZZ T1;

   xdouble *mu_k = mu[k];

   if (st < k) {
      for (i = 1; i < st; i++)
         buf[i] = mu_k[i]*c[i];
   }

   for (j = st; j <= k-1; j++) {
      if (b[k]*b[j] < NTL_FDOUBLE_PRECISION*NTL_FDOUBLE_PRECISION) {
         double z = 0;
         xdouble *B1_k = B1[k];
         xdouble *B1_j = B1[j];

         for (i = 1; i <= n; i++)
            z += B1_k[i].x * B1_j[i].x;

         s = z;
      }
      else {
         s = InnerProduct(B1[k], B1[j], n);
   
         if (s*s <= b[k]*b[j]/bound) {
            InnerProduct(T1, B(k), B(j));
            conv(s, T1);
         }
      }

      xdouble *mu_j = mu[j];

      t1 = 0;
      for (i = 1; i <= j-1; i++)
         MulAdd(t1, t1, mu_j[i], buf[i]);

      mu_k[j] = (buf[j] = (s - t1))/c[j];
   }

   s = 0;
   for (j = 1; j <= k-1; j++)
      MulAdd(s, s, mu_k[j], buf[j]);

   c[k] = b[k] - s;
}

double Deriv1 (edge *TheEdge, position *ThePosition, line *TheLine,
	       int index)
{
  double v[3], newv[3], dv[3], hat[3];

  if (!TheEdge || !ThePosition || !TheLine) return (0.0);

  switch (index) {
    case THETA :
      hat[0] = 0.0;  hat[1] = -1.0;  hat[2] = 0.0;
      RotateQuaternion (dv, TheLine->q, hat);
      RotateQuaternion (newv, ThePosition->q, dv);
      return (InnerProduct (newv, TheEdge->m));
      break;

    case OMEGA :
      hat[0] = 1.0;  hat[1] = 0.0;  hat[2] = 0.0;
      RotateQuaternion (dv, TheLine->q, hat);
      RotateQuaternion (newv, ThePosition->q, dv);
      return (InnerProduct (newv, TheEdge->m));
      break;

    case ALPHA :
      ComputeV (v, *TheLine);
      RotateQuaternion (newv, ThePosition->q, v);
      DerivR (dv, newv, 'x');
      return (InnerProduct (dv, TheEdge->m));
      break;

    case BETA :
      ComputeV (v, *TheLine);
      RotateQuaternion (newv, ThePosition->q, v);
      DerivR (dv, newv, 'y');
      return (InnerProduct (dv, TheEdge->m));
      break;

    case GAMMA :
      ComputeV (v, *TheLine);
      RotateQuaternion (newv, ThePosition->q, v);
      DerivR (dv, newv, 'z');
      return (InnerProduct (dv, TheEdge->m));
      break;

    default :
      pexit ("Bad parameters to Deriv1");
      break;
  }

  return 0.0;
}

double SBVAR::LogLikelihood(void)
{
  double log_likelihood=log_likelihood_constant + lambda_T*LogAbsDeterminant(A0);
  TDenseVector a0(n_vars), aplus(n_predetermined);

  for (int i=n_vars-1; i >= 0; i--)
    {
      a0.RowVector(A0,i);
      aplus.RowVector(Aplus,i);
      log_likelihood+=-0.5*(InnerProduct(a0,a0,YY) - 2.0*InnerProduct(aplus,a0,XY) + InnerProduct(aplus,aplus,XX));
    }

  return log_likelihood;
}

/*
   Computes the log of the integral of the tempered posterior.

         Integrate( (p(Y|Theta)*p(Theta))^(1/K) dTheta )
*/
double SBVAR_symmetric_linear::LogPosteriorIntegral(TDenseVector p, int ndraws, int thin, int burn_in)
{
  if (ndraws <= 0) throw dw_exception("PosteriorIntegral(): number of draws must be postive");
  if (thin <= 0) throw dw_exception("PosteriorIntegral(): thinning factor must be positive");
  if (p.dim != NumberParameters()) throw dw_exception("PosteriorIntegral(): Incorrect number of parameters");

  if (!simulation_info_set) SetSimulationInfo();

  SetParameters(p.vector);

  double integral=log_likelihood_constant + log_prior_constant;
  for (int i=0; i < n_vars; i++)
    integral+=0.5*dim_g[i]*1.837877066409345 + LogAbsDeterminant(Simulate_SqrtH[i]);  // 1.837877066409345 = log(2*pi)

  integral+=lambda_T*LogAbsDeterminant(A0);
  for (int i=0; i < n_vars; i++)
    {
      TDenseVector x=InverseMultiply(Simulate_SqrtS[i],p.SubVector(begin_b[i],begin_b[i]+dim_b[i]-1));
      integral+=-0.5*lambda_T*InnerProduct(x,x);
    }

  for (int i=0; i < n_vars; i++)
    integral-=LogConditionalA0_gibbs(p,i,ndraws,thin,burn_in);

  return integral;
}

    wrapper(size_t n, params prm = params(),
            const backend_params &bprm = backend_params(),
            const InnerProduct &inner_product = InnerProduct()
            )
        : s(prm.get("type", runtime::solver::bicgstab)), handle(0)
    {
        if (!prm.erase("type")) AMGCL_PARAM_MISSING("type");

        switch(s) {

#define AMGCL_RUNTIME_SOLVER(type) \
            case type: \
                handle = static_cast<void*>(new amgcl::solver::type<Backend, InnerProduct>(n, prm, bprm, inner_product)); \
                break

            AMGCL_RUNTIME_SOLVER(cg);
            AMGCL_RUNTIME_SOLVER(bicgstab);
            AMGCL_RUNTIME_SOLVER(bicgstabl);
            AMGCL_RUNTIME_SOLVER(gmres);
            AMGCL_RUNTIME_SOLVER(lgmres);
            AMGCL_RUNTIME_SOLVER(fgmres);
            AMGCL_RUNTIME_SOLVER(idrs);

#undef AMGCL_RUNTIME_SOLVER

            default:
                throw std::invalid_argument("Unsupported solver type");
        }
    }

    double GillespiePetzold_Stepsize(Vector& x ,
                          Vector& a ,
                          double & a0 ,
                          const Matrix& nu ,
                          double tau,
                          double eps, 
                          PropensityJacobianFunc jacobian)
    {
	static Matrix f(x.Size(), x.Size()); 
	static Vector mu = a; 
	static Vector sigma = a; 
	static double r1, r2; 
	static double tau1, tau2; 

	f = jacobian(x)*nu; 
	mu = f*a; 
	sigma = InnerProduct(f, f)*a; 
	r1 = mu.Norm(); //Max norm
	r2 = sigma.Norm(); // Max norm
	tau1 = eps * a0; 
	tau2 = tau1 * tau1; 
	if (r1 > 0) 
	  tau1 = tau1/r1; 
	else 
	  tau1 = 100.; 
	if (r2 > 0) 
	  tau2 = tau2/r2; 
	else 
	  tau2 = 100.; 
 
	if (tau1 > tau2) 
	  return tau2; 
	else 
	  return tau1; 
    }

NTL_START_IMPL

static
long CharPolyBound(const mat_ZZ& a)
// This bound is computed via interpolation
// through complex roots of unity.

{
   long n = a.NumRows();
   long i;
   ZZ res, t1, t2;

   set(res);

   for (i = 0; i < n; i++) {
      InnerProduct(t1, a[i], a[i]);
      abs(t2, a[i][i]);
      mul(t2, t2, 2);
      add(t2, t2, 1);
      add(t1, t1, t2);
      if (t1 > 1) {
         SqrRoot(t1, t1);
         add(t1, t1, 1);
      }
      mul(res, res, t1);
   }

   return NumBits(res);
}

TInt DoCP0Test(TAny* aPtr)
	{
	TUint seed[2];
	seed[0]=(TUint)aPtr;
	seed[1]=0;
	TInt16 vec1[128];
	TInt16 vec2[128];
	TInt run;
	for (run=0; run<100000; ++run)
		{
		TInt n=(Random(seed)&63)+64;		// vector length
		TInt i;
		for (i=0; i<n; ++i)
			{
			vec1[i]=(TInt16)(Random(seed)&0xffff);
			vec2[i]=(TInt16)(Random(seed)&0xffff);
			TInt64 result, result2;
			InnerProduct(result,vec1,vec2,n);
			InnerProduct2(result2,vec1,vec2,n);
			if (result != result2)
				{
				User::Panic(_L("ERROR"),run);
				}
			}
		}
	return 0;
	}

 double ScalarFiniteElement<D> :: 
 Evaluate (const IntegrationPoint & ip, FlatVector<double> x) const
 {
   VectorMem<20, double> shape(ndof);
   CalcShape (ip, shape);
   return InnerProduct (shape, x);
 }  

void ProjectPowers(vec_zz_p& x, const vec_zz_p& a, long k,
                   const zz_pXArgument& H, const zz_pXModulus& F)

{
   long n = F.n;

   if (a.length() > n || k < 0 || NTL_OVERFLOW(k, 1, 0))
      Error("ProjectPowers: bad args");

   long m = H.H.length()-1;
   long l = (k+m-1)/m - 1;

   zz_pXMultiplier M;
   build(M, H.H[m], F);

   vec_zz_p s(INIT_SIZE, n);
   s = a;
   StripZeroes(s);

   x.SetLength(k);

   for (long i = 0; i <= l; i++) {
      long m1 = min(m, k-i*m);
      zz_p* w = &x[i*m];
      for (long j = 0; j < m1; j++)
         InnerProduct(w[j], H.H[j].rep, s);
      if (i < l)
         UpdateMap(s, s, M, F);
   }
}