bool NonlinearDriver::solutionNorms (const TimeDomain& time,
                                     double zero_tol, std::streamsize outPrec)
{
  if (msgLevel < 0 || solution.empty()) return true;

  const size_t nsd = model.getNoSpaceDim();

  size_t iMax[nsd];
  double dMax[nsd];
  double normL2 = model.solutionNorms(solution.front(),dMax,iMax);

  RealArray RF;
  bool haveReac = model.getCurrentReactions(RF,solution.front());

  Vectors gNorm;
  if (calcEn)
  {
    model.setMode(SIM::RECOVERY);
    model.setQuadratureRule(opt.nGauss[1]);
    if (!model.solutionNorms(time,solution,gNorm))
      gNorm.clear();
  }

  if (myPid > 0) return true;

  std::streamsize stdPrec = outPrec > 0 ? IFEM::cout.precision(outPrec) : 0;
  double old_tol = utl::zero_print_tol;
  utl::zero_print_tol = zero_tol;

  IFEM::cout <<"  Primary solution summary: L2-norm            : "
             << utl::trunc(normL2);

  for (unsigned char d = 0; d < nsd; d++)
    if (utl::trunc(dMax[d]) != 0.0)
      IFEM::cout <<"\n                            Max "<< char('X'+d)
                 <<"-displacement : "<< dMax[d] <<" node "<< iMax[d];

  if (haveReac)
  {
    IFEM::cout <<"\n  Total reaction forces: Sum(R) =";
    for (size_t i = 1; i < RF.size(); i++)
      IFEM::cout <<" "<< utl::trunc(RF[i]);
    if (utl::trunc(RF.front()) != 0.0)
      IFEM::cout <<"\n  displacement*reactions: (R,u) = "<< RF.front();
  }

  if (!gNorm.empty())
    this->printNorms(gNorm.front(),IFEM::cout);

  IFEM::cout << std::endl;
  utl::zero_print_tol = old_tol;
  if (stdPrec > 0) IFEM::cout.precision(stdPrec);
  return true;
}

  static boost::python::list
  get_identifiers(const Vectors& vec)
  {
    boost::python::list l;

    int nb_vec = vec.get_nb_vector();
    for (int v = 0; v < nb_vec; v++)
      {
        l.append(vec.get_identifier(v));
      }

    return l;
  }

 // check
 static bool
 check(Vectors& v)
 {
   StatError error;
   bool ret = v.check(error);
   return ret;
 }

  virtual bool assembleSystem(const TimeDomain& time,
                              const Vectors& prevSol,
                              bool newLHSmatrix, bool)
  {
    const double M = 10.0;   // Mass of the oscillator
    const double K = 1000.0; // Stiffness of the oscillator
    const double F = 1.0;    // External load (constant)

    myEqSys->initialize(newLHSmatrix);

    bool ok;
    if (myProblem->getMode() == SIM::MASS_ONLY) {
      ElmMats elm;
      elm.resize(1,1); elm.redim(1);
      elm.A[0].fill(M); // Mass Matrix
      ok = myEqSys->assemble(&elm,1);
    }
    else {
      const double* intPrm = static_cast<Problem*>(myProblem)->getIntPrm();
      NewmarkMats elm(intPrm[0],intPrm[1],intPrm[2],intPrm[3]);
      elm.resize(3,1); elm.redim(1); elm.vec.resize(3);
      elm.setStepSize(time.dt,0);
      elm.A[1].fill(M); // Mass matrix
      elm.A[2].fill(K); // Stiffness matrix
      elm.b[0] = -K*prevSol.front(); // Elastic forces
      for (int i = 0; i < 3; i++) elm.vec[i] = prevSol[i];
      ok = myEqSys->assemble(&elm,1);
    }

    // Add in the external load
    ok &= mySam->assembleSystem(*myEqSys->getVector(),&F,1);

    return ok && myEqSys->finalize(newLHSmatrix);
  }

int main(){
	Vectors<int>* newVector = new Vectors<int>(10,10,0);
	for (int i=0; i<10; ++i){
		for (int j=0; j<10; ++j){
			newVector->insert(i, j, 5);
		}
	}
	newVector->print();
	
	for (int i=0; i<10; ++i){
		for (int j=0; j<10; ++j){
			newVector->remove(i, j);
		}
	}
	newVector->print();
	delete newVector;
	return 0;
}

double SIMLinElKL::externalEnergy (const Vectors& psol) const
{
  double energy = this->SIMbase::externalEnergy(psol);

  // External energy from the nodal point loads
  for (size_t i = 0; i < myLoads.size(); i++)
    energy += myLoads[i].pload * psol.front()(myLoads[i].inod);

  return energy;
}

TEST(Vectors, Insert) {
	Vectors<int>* newInt = new Vectors<int>(10, 10, 0);
	newInt->insert(4,5,5);
	EXPECT_EQ(5, newInt->access(4,5));
	delete newInt;

	Vectors<std::string>* newString = new Vectors<std::string>(9, 18, "empty");
	newString->insert(4,5,"hello");
	EXPECT_EQ("hello", newString->access(4,5));
	delete newString;

	Vectors<double>* newDouble = new Vectors<double>(6, 7, 0.0);
	newDouble->insert(4,5,5.5);
	EXPECT_EQ(5.5, newDouble->access(4,5));
	delete newDouble;
}

  static bool
  select_step(Vectors &input, int variable, double step)
  {
    StatError error;
    bool ret;

    ret = input.select_step(error, variable, step);
    if (!ret)
      stat_tool::wrap_util::throw_error(error);
    return ret;
  }

  static void
  file_ascii_data_write(const Vectors& d, const char* path, bool exhaustive)
  {
    bool result = true;
    StatError error;

    result = d.ascii_data_write(error, path, exhaustive);
    if (!result)
      stat_tool::wrap_util::throw_error(error);

  }

  static std::string
  ascii_data_write(const Vectors& d, bool exhaustive)
  {
    std::stringstream s;
    std::string res;

    d.ascii_data_write(s, exhaustive);
    res = s.str();

    return res;

  }

 static string
 rank_correlation_computation(const Vectors& input, int icorrel_type, const string &filename)
 {
   StatError error;
   std::stringstream os;
   bool ret;
   correlation_type correl_type = correlation_type(icorrel_type);
   
   ret = input.rank_correlation_computation(error, os, correl_type, filename.c_str());
   //std::cout << os.str()<<endl;
   return os.str();
 }

double SIMKLShell::externalEnergy (const Vectors& u,
                                   const TimeDomain& time) const
{
  double energy = this->SIMbase::externalEnergy(u,time);

  // External energy from the nodal point loads
  const int* madof = mySam->getMADOF();
  for (const PointLoad& load : myLoads)
    if (load.ldof.second > 0)
    {
      int idof = madof[load.ldof.first-1] + load.ldof.second-1;
      energy += (*load.p)(time.t) * u.front()(idof);
    }
    else if (load.ldof.second < 0) // This is an element point load
    {
      Vector v = this->SIMgeneric::getSolution(u.front(),load.xi,0,load.patch);
      if (-load.ldof.second <= (int)v.size())
        energy += (*load.p)(time.t) * v(-load.ldof.second);
    }

  return energy;
}

  static string
  contingency_table(const Vectors& v, int variable1, int variable2,
      const string& filename, int iformat)
  {
    StatError error;
    std::stringstream s;
    bool ret;
    output_format format = output_format(iformat);   
    
    ret = v.contingency_table(error, s, variable1, variable2, filename.c_str(),
        format);

    if (!ret)
      stat_tool::wrap_util::throw_error(error);

    return s.str();
  }

  static string
  variance_analysis(const Vectors& v, int class_variable,
      int response_variable, int response_type, const string& filename,
      int iformat)
  {
    StatError error;
    std::stringstream s;
    bool ret;
    output_format format = output_format(iformat);

    ret = v.variance_analysis(error, s, class_variable, response_variable,
        response_type, filename.c_str(), format);

    if (!ret)
      stat_tool::wrap_util::throw_error(error);

    return s.str();
  }

bool ASMunstruct::refine (const LR::RefineData& prm,
                          Vectors& sol, const char* fName)
{
  PROFILE2("ASMunstruct::refine()");

  if (!geo)
    return false;
  else if (shareFE && !prm.refShare)
  {
    nnod = geo->nBasisFunctions();
    return true;
  }

  // to pick up if LR splines get stuck while doing refinement,
  // print entry and exit point of this function

  double beta         = prm.options.size() > 0 ? prm.options[0]/100.0 : 0.10;
  int    multiplicity = prm.options.size() > 1 ? prm.options[1]       : 1;
  if (multiplicity > 1)
    for (int d = 0; d < geo->nVariate(); d++) {
      int p = geo->order(d) - 1;
      if (multiplicity > p) multiplicity = p;
    }

  enum refinementStrategy strat = LR_FULLSPAN;
  if (prm.options.size() > 2)
    switch (prm.options[2]) {
    case 1: strat = LR_MINSPAN; break;
    case 2: strat = LR_STRUCTURED_MESH; break;
    }

  bool linIndepTest   = prm.options.size() > 3 ? prm.options[3] != 0 : false;
  int  maxTjoints     = prm.options.size() > 4 ? prm.options[4]      : -1;
  int  maxAspectRatio = prm.options.size() > 5 ? prm.options[5]      : -1;
  bool closeGaps      = prm.options.size() > 6 ? prm.options[6] != 0 : false;

  char doRefine = 0;
  if (!prm.errors.empty())
    doRefine = 'E'; // Refine based on error indicators
  else if (!prm.elements.empty())
    doRefine = 'I'; // Refine the specified elements

  if (doRefine) {

    std::vector<int> nf(sol.size());
    for (size_t j = 0; j < sol.size(); j++)
      if (!(nf[j] = LR::extendControlPoints(geo,sol[j],this->getNoFields(1))))
        return false;

    // set refinement parameters
    if (maxTjoints > 0)
      geo->setMaxTjoints(maxTjoints);
    if (maxAspectRatio > 0)
      geo->setMaxAspectRatio((double)maxAspectRatio);
    geo->setCloseGaps(closeGaps);
    geo->setRefMultiplicity(multiplicity);
    geo->setRefStrat(strat);

    // do actual refinement
    if (doRefine == 'E')
      geo->refineByDimensionIncrease(prm.errors,beta);
    else if (strat == LR_STRUCTURED_MESH)
      geo->refineBasisFunction(prm.elements);
    else
      geo->refineElement(prm.elements);

    geo->generateIDs();
    nnod = geo->nBasisFunctions();

    for (int i = sol.size()-1; i >= 0; i--) {
      sol[i].resize(nf[i]*geo->nBasisFunctions());
      LR::contractControlPoints(geo,sol[i],nf[i]);
    }
  }

  if (fName)
  {
    char fullFileName[256];

    strcpy(fullFileName, "lrspline_");
    strcat(fullFileName, fName);
    std::ofstream lrOut(fullFileName);
    lrOut << *geo;
    lrOut.close();

    LR::LRSplineSurface* lr = dynamic_cast<LR::LRSplineSurface*>(geo);
    if (lr) {
      // open files for writing
      strcpy(fullFileName, "param_");
      strcat(fullFileName, fName);
      std::ofstream paramMeshFile(fullFileName);

      strcpy(fullFileName, "physical_");
      strcat(fullFileName, fName);
      std::ofstream physicalMeshFile(fullFileName);

      strcpy(fullFileName, "param_dot_");
      strcat(fullFileName, fName);
      std::ofstream paramDotMeshFile(fullFileName);

//.........这里部分代码省略.........

//.........这里部分代码省略.........
        if(!mxIsClass(tmp,"int32")) {
            mexErrMsgTxt("input.classname must be int32\n");
        }
        else    {
            int *pVECTOR_CLASS = (int*)mxGetData(tmp);
            VECTOR_CLASS = (VectorClassNames)(*pVECTOR_CLASS);
        }
    }
    
    //mexPrintf("DISTANCE_FCN=%d\n",DISTANCE_FCN);
    
    /* Get dimensions of second field of input */
    
    tmp=prhs[1];
    mwSize ndims_in=mxGetNumberOfDimensions(tmp);
    const mwSize* dims_in=mxGetDimensions(tmp);
    if(!mxIsClass(tmp,mxIsClassName)) {
        mexErrMsgTxt("Second input must be double\n");
    }
    
    
//    mexPrintf("ndims_in=%d\n",ndims_in);
    
    int N=dims_in[ndims_in-1];
    dim=1;
    for(int i=0;i<ndims_in-1;i++) {
        dim*=dims_in[i];
    }
    
//    mexPrintf("dim=%d N=%d\n",dim,N);
    
    // Create Vectors
    REAL* X=(REAL*)mxGetData(tmp);
    Vectors *vectors;
    switch (VECTOR_CLASS)   {
        case VECTORS:
            vectors = new Vectors(X,N,dim,DISTANCE_FCN);
            break;
        default:
            mexErrMsgTxt("\n Invalid classname for Vectors.");
            break;
    }
    
    if( N==1 ) {
        NTHREADS = 0;
    }
    ThreadsWithCounter threads(NTHREADS);
    
    // Construct covertree
    SegList<DLPtrListNode<CoverNode> > seglist(1024);
    const Vector* vector=(vectors->next());
    Cover cover(vector,seglist,numlevels,minlevel);
    EnlargeData enlargedata(&threads,BLOCKSIZE,vectors->getRemaining());
//    Timer timer;
//    timer.on();
    cover.enlargeBy(enlargedata,*vectors);
//    timer.off();
//    timer.printOn(cout);
    
    
    /* Create matrix for the return argument. */
    plhs[0] = mxCreateStructMatrix(1, 1, 5, fnames_out);
    
    mxArray* fout;
    
    dims[0]=1;

int main(void) {

  int var, i, j, nb_variable, nb_component;
  const int nb_vector = 500;
  double *pweight = NULL;
  bool *fparam = NULL;
  int* perm;
  StatError error;
  const char * mixpath= "./tmp.mix", * spmixpath= "./tmp_sp.mix";
  const char * gnupath = "./tmp_mix", * gnu_datapath = "./tmp_mix_data";
  const char * margpath= "./marg_mix", * gnu_tmppath = "./tmp_mix_d";
  const char * np_modelpath= "./np_model.mix", * gnu_tmpnppath = "./tmp_mix_d";
  const char * entropy_data= "cluster_vectors.vec";
  const char * output_path = "tmp_entropy_output.vec";
  Distribution *marginal = NULL;
  DiscreteDistributionData *marginal_histo = NULL;
  Vectors *vec = NULL;
  MultivariateMixture *mv1 = NULL, *mv_cp = NULL;
  MultivariateMixture *mv_np1 = NULL, *mv_np_estim = NULL;
  MultivariateMixture *mv_estim = NULL;
  MultivariateMixtureData *mv_data = NULL, *cluster = NULL;
  DiscreteParametric **dt1 = NULL, **dt2 = NULL;
  DiscreteParametricProcess **ppcomponent = NULL;


  // constructors of Mv_Mixture
  mv1 = new MultivariateMixture();

  // destructor of MultivariateMixture
  delete mv1;
  mv1= NULL;

  nb_variable = 2;
  nb_component = 3;

  dt1 = new DiscreteParametric*[nb_component];
  dt2 = new DiscreteParametric*[nb_component];
  pweight = new double[nb_component];
  ppcomponent = new DiscreteParametricProcess*[nb_variable];

  pweight[0] = 0.1;
  pweight[1] = 0.2;
  pweight[2] = 0.7;



  dt1[0] = new DiscreteParametric(0, BINOMIAL, 2, 12, D_DEFAULT, 0.1);
  dt1[1] = new DiscreteParametric(0, BINOMIAL, 0, 10, D_DEFAULT, 0.5);
  dt1[2] = new DiscreteParametric(0, BINOMIAL, 3, 10, D_DEFAULT, 0.8);

  dt2[0] = new DiscreteParametric(0, POISSON, 2, I_DEFAULT, 8.0, D_DEFAULT);
  dt2[1] = new DiscreteParametric(0, POISSON, 4, I_DEFAULT, 5.0, D_DEFAULT);
  dt2[2] = new DiscreteParametric(0, POISSON, 0, I_DEFAULT, 2.0, D_DEFAULT);

  cout << "Observation distributions for variable 1:" << endl;
  for (i = 0; i < nb_component; i++) {
    dt1[i]-> ascii_print(cout);
  }

  ppcomponent[0] = new DiscreteParametricProcess(nb_component, dt1);
  ppcomponent[1] = new DiscreteParametricProcess(nb_component, dt2);

  for (i = 0; i < nb_component; i++) {
    delete dt1[i];
    dt1[i] = NULL;
    delete dt2[i];
    dt2[i] = NULL;
  }

  cout << endl;

  mv1 = new MultivariateMixture(nb_component, pweight, nb_variable, ppcomponent, NULL);

  cout << "Mixture of " << nb_component << " components with " <<
    nb_variable << " variables:" << endl;

  mv1->ascii_write(cout, true);
  cout << endl;

  // copy
  mv_cp = new MultivariateMixture(*mv1);
  cout << "Copy constructor of MultivariateMixture: " << endl;
  mv_cp->ascii_write(cout);
  cout << endl;

  // destructor of MultivariateMixture
  delete mv_cp;
  mv_cp= NULL;

  delete mv1;
  mv1= NULL;

  cout << "MultivariateMixture_building (print into file " << mixpath << "): " << endl;

  mv1 = multivariate_mixture_building(error , nb_component ,
                                      nb_variable, pweight, ppcomponent, NULL);

  if (mv1 == NULL)
    cout << error;
  else {
//.........这里部分代码省略.........

void __stdcall ClearVectors()
{
  for (Vectors::size_type i = 0; i < vectors.size(); ++i)
    delete[] vectors[i];
  vectors.clear();
}

 static void plot_write(const Vectors &input, const std::string& prefix,
 const std::string& title)
 {
     StatError error;
     input.plot_write(error, prefix.c_str(), title.c_str());
 }

bool Elasticity::evalSol (Vector& s, const Vectors& eV, const FiniteElement& fe,
                          const Vec3& X, bool toLocal, Vec3* pdir) const
{
  if (eV.empty())
  {
    std::cerr <<" *** Elasticity::evalSol: No solutions vector."<< std::endl;
    return false;
  }
  else if (!eV.front().empty() && eV.front().size() != fe.dNdX.rows()*nsd)
  {
    std::cerr <<" *** Elasticity::evalSol: Invalid displacement vector."
	      <<"\n     size(eV) = "<< eV.front().size() <<"   size(dNdX) = "
	      << fe.dNdX.rows() <<","<< fe.dNdX.cols() << std::endl;
    return false;
  }

  // Evaluate the deformation gradient, dUdX, and/or the strain tensor, eps
  Matrix Bmat;
  Tensor dUdX(nDF);
  SymmTensor eps(nsd,axiSymmetry);
  if (!this->kinematics(eV.front(),fe.N,fe.dNdX,X.x,Bmat,dUdX,eps))
    return false;

  // Add strains due to temperature expansion, if any
  double epsT = this->getThermalStrain(eV.back(),fe.N,X);
  if (epsT != 0.0) eps -= epsT;

  // Calculate the stress tensor through the constitutive relation
  Matrix Cmat;
  SymmTensor sigma(nsd, axiSymmetry || material->isPlaneStrain()); double U;
  if (!material->evaluate(Cmat,sigma,U,fe,X,dUdX,eps))
    return false;
  else if (epsT != 0.0 && nsd == 2 && material->isPlaneStrain())
    sigma(3,3) -= material->getStiffness(X)*epsT;

  Vec3 p;
  bool havePval = false;
  if (toLocal && wantPrincipalStress)
  {
    // Calculate principal stresses and associated direction vectors
    if (sigma.size() == 4)
    {
      SymmTensor tmp(2); tmp = sigma; // discard the sigma_zz component
      havePval = pdir ? tmp.principal(p,pdir,2) : tmp.principal(p);
    }
    else
      havePval = pdir ? sigma.principal(p,pdir,2) : sigma.principal(p);

    // Congruence transformation to local coordinate system at current point
    if (locSys) sigma.transform(locSys->getTmat(X));
  }

  s = sigma;

  if (toLocal)
    s.push_back(sigma.vonMises());

  if (havePval)
  {
    s.push_back(p.x);
    s.push_back(p.y);
    if (sigma.dim() == 3)
      s.push_back(p.z);
  }

  return true;
}