/*!
 * Initializes the rhs_CNLS class  
 */
void rhs_CNLS::postprocess(input& dat){
    rhs::postprocess(dat);
    NUM_TIME_STEPS = dimension/2;
    if(NUM_TIME_STEPS*2 != dimension){
        err("dimension not even, which is required for rhs_CNLS", 
                "rhs_CNLS::postprocess", "rhs/rhs_CNLS.cpp", FATAL_ERROR);
    }
    dat.retrieve(LENGTH_T, "t_int", this);
    if(LENGTH_T <= 0){
        std::string errmess = "t_int is invalid, must be >= 0";
        err(errmess, "rhs_CNLS::postprocess", "rhs/rhs_CNLS.cpp",
                dat["t_int"], FATAL_ERROR);
    }
    dt = LENGTH_T/NUM_TIME_STEPS;
    dat.retrieve(g0, "g0", this);
    dat.retrieve(e0, "e0", this);
    memp.create(NUM_TIME_STEPS, &u1, &u2, &comp_in, &comp_in_r, &comp_out, &comp_out_r, &sq1, &sq2, &k, &ksq);
    //create k values
    double mulval=(2.0*PI/LENGTH_T)*(NUM_TIME_STEPS/2.0);
    for(size_t i=0; i<NUM_TIME_STEPS/2; i++){
        k[i] = mulval * (2.0*i/(1.0*NUM_TIME_STEPS));
        ksq[i] = k[i]*k[i];
    }
    for(size_t i=NUM_TIME_STEPS/2; i<NUM_TIME_STEPS; i++){
        k[i] = mulval * 2.0*((int)i-(int)NUM_TIME_STEPS)/(NUM_TIME_STEPS*1.0);
        ksq[i] = k[i]*k[i];
    }

}

template<typename String, typename Iter> inline bool _parse_codepoint(String& out, input<Iter>& in) {
  int uni_ch;
  if ((uni_ch = _parse_quadhex(in)) == -1) {
    return false;
  }
  if (0xd800 <= uni_ch && uni_ch <= 0xdfff) {
    if (0xdc00 <= uni_ch) {
  // a second 16-bit of a surrogate pair appeared
  return false;
    }
    // first 16-bit of surrogate pair, get the next one
    if (in.getc() != '\\' || in.getc() != 'u') {
  in.ungetc();
  return false;
    }
    int second = _parse_quadhex(in);
    if (! (0xdc00 <= second && second <= 0xdfff)) {
  return false;
    }
    uni_ch = ((uni_ch - 0xd800) << 10) | ((second - 0xdc00) & 0x3ff);
    uni_ch += 0x10000;
  }
  if (uni_ch < 0x80) {
    out.push_back(uni_ch);
  } else {
    if (uni_ch < 0x800) {
  out.push_back(0xc0 | (uni_ch >> 6));
    } else {
  if (uni_ch < 0x10000) {

/*!
 * Performs the postprocessing for toroidal,
 * which is setting the number of iterations that are to be performed
 * @param dat Map containing the input values
 * \sa controller::postprocess, item_dim::postprocess
 */
void toroidal::postprocess(input& dat){
    controller::postprocess(dat);
    num_int=0;
    int _iterations;



    dat.retrieve(_iterations, "iterations", this);
    dat.retrieve(initial_inc, "initial_inc", this);
    dat.retrieve(mul_fac, "mul_fac", this);


    iterations = _iterations;
    if(mul_fac==0){
        err("Multiply factor, mul_fac, must not be equal to zero", 
               "toroidal::postprocess", "controller/toroidal.cpp",
               FATAL_ERROR);
    }
    if(initial_inc==0){
        err("The initial increment, initial_inc, must not be equal to zero", 
               "toroidal::postprocess", "controller/toroidal.cpp",
               FATAL_ERROR);
    }
    //temporary hack to help with c_elegans
    holder->index = index*iterations;
    //find the controllers place in the number of iterations
}

/*!
 * This function process data that has been generated from an input file
 * All of the names returned in the dependencies file exist and have been processed
 * @param dat The data structure containing the variable names
 */
void example_rhs::postprocess(input& dat){
    //warning so nobody ever actually constructs one of these
    //please don't have warnings that always throw in your code ever
    err("Example_rhs created, mysteriously fails on some architectures seemingly from AVX instructions. Also does nothing useful", "rhs::create", "rhs/rhs.cpp", WARNING);

    rhs::postprocess(dat);//always postprocess the parent class first
    //Any class that inherits from rhs, or item_dim, has access to a variable
    //called dimension, which represents the dimension of the problem at hand
    //dimension is initialized in the item_dep, which is why the parent postprocessing
    //must be called first
    

    //There are two types of errors, warning and fatal errors.
    //Warnings print something to the screen to warn the user,
    //but do not stop the program. A fatal error causes an immediate exit
    //These are caused by the function err. Details can be found in the documentation
    //
    //Lets throw a fatal error if the dimension is less than 10
    if(dimension > 10){
        err("dimension > 10, which is required for example_rhs", "example_rhs::postprocess",
                "rhs/example_rhs.cpp", FATAL_ERROR);
    }
    
    //Lets also throw an error if the dimension is equal to 5
    if(dimension==5){
        err("dimension=5, which gets a warning from example_rhs", "example_rhs::postprocess",
                "rhs/example_rhs.cpp", WARNING);
    }
    //The map contains item*, which point to an item of arbitrary type.
    //I plan to implement type checking of the dependencies.
    //But sometime later
    //
    //To retrieve a value, you can call the retrieve function of the variable at hand
    //you just pass the name of parameter that is being retrieved to the map, which returns
    //an item*. Then this pointer is called to retrieve the value, and is passed the address
    //of val1

    dat.retrieve(val1, "val1", this);
    dat.retrieve(val2, "val2", this);
    dat.retrieve(random_info, "random_info", this);

    //now, we are going to allocate some memory to something
    //This may be useful for storing temporary calculations during the RHS.
    //Inheriting from item_dim provides access to a memory pool, memp.
    //just pass memp.create the dimension of the problem and the
    //addresses of the pointers and the alignment, allocation, and deallocation is
    //taken care of!
    //
    //memp.create(dimension, &ptr1, &ptr2, &ptr3, etc)
    //for custom alignment(standard is 32 byte)
    //memp.create(alignment, dimension, &ptr1, &prt2, etc)
    //
    memp.create(32, dimension, &value_holder);
    for(size_t i = 0; i < dimension; i++){
        value_holder[i] = Id*(double)i*val1 + (dimension-i)*val2;
    }
}

inline void read_input(input &inp, Problems &p)
{
  unsigned short u, v;

  p.reset();

  inp.read_short(u);

  for(unsigned short i = 0; i < u; ++i)
  {
    inp.read_short(v);
    p.set(v);
  }
}

int copy_file_thru_input::inget()
{
  if (!in)
    return EOF;
  else
    return in->get();
}

	output fc_rnn::execute(input const& in)
	{
		// Set activation of input neurons
		auto const num_input = in.size();
		for(size_t n = 0; n < num_input; ++n)
		{
			vInput[n] = in[n];
		}

		// Summation for hidden neurons
		Eigen::VectorXd vHiddenSums =
			wmInput * vInput +
			wmHidden * vHidden;
		// Transfer function
		vHidden =
			evaluate(af_hidden, vHiddenSums.array());

		// TODO: Maybe should just store as a single vector?
		Eigen::VectorXd joined(input_layer_count() + hidden_count());
		joined << vInput, vHidden;
		Eigen::VectorXd vOutputSums =
			wmOutput * joined;
		Eigen::VectorXd vOutput =
			evaluate(af_output, vOutputSums.array());

		// Return the output values
		output out{ output_count() };
		std::copy(vOutput.data(), vOutput.data() + output_count(), out.begin());
		return out;
	}

/**
 * Calculate the nuclear-nuclear potential energy of a problem
 * @param problem the problem to be solved
 */
double MxElem::NuclPotEn(input & problem){

   double energy = 0;
   int Ncores = problem.gNcores();

   for (int i=0; i<Ncores; i++){

      R Ri(*(problem.gvector(i)));
      int Zi = problem.gcore(i);

      for (int j=i+1; j<Ncores; j++){

         energy += Zi * problem.gcore(j) / sqrt(Ri.DistanceSquared(*(problem.gvector(j))));
      }
   }

   return energy;

}

void example_integrator_tmpl<T>::postprocess(input& in){
    //perform postprocessing of integrator class
    //note how we skip example_integrator in this chain since example integrator is
    //only a proxy that performs type erasure
    integrator::postprocess(in);
    item* some_class;
    in.retrieve(some_class, "test_class", this);
    //this allows us to have a consistent representation in the input file.
    //If the variable in question is not part of the hardcore numerical analysis code,
    //it may be better to just use a double and not deal with this in the postprocessing
    double _rval1 = 0; 
    in.retrieve(_rval1, "rval1", this);
    rval1=_rval1;

    //note how we pass a default parameter since rval2 isn't necesarily going to exist
    in.retrieve(rval2, "rval2", this, 0);
    in.retrieve(unsigned_var, "unsigned_var", this);
    in.retrieve(something, "something", this);
    update();
}

 template<typename Iter> inline int _parse_quadhex(input<Iter> &in) {
   int uni_ch = 0, hex;
   for (int i = 0; i < 4; i++) {
     if ((hex = in.getc()) == -1) {
   return -1;
     }
     if ('0' <= hex && hex <= '9') {
   hex -= '0';
     } else if ('A' <= hex && hex <= 'F') {
   hex -= 'A' - 0xa;
     } else if ('a' <= hex && hex <= 'f') {
   hex -= 'a' - 0xa;
     } else {
   in.ungetc();
   return -1;
     }
     uni_ch = uni_ch * 16 + hex;
   }
   return uni_ch;
 }

void n_pulse_score::postprocess(input& invals){
    objective::postprocess(invals);
    int _n_pulse = 0;
    invals.retrieve(_n_pulse, "num_pulses", this);
    n_pulse = _n_pulse;
    if(dimension%n_pulse){
        err("n_pulse_score requires a dimension divisible by the number of pulses", 
                "bi_pulse_score::postprocess", "objective/bi_pulse_score.cpp", FATAL_ERROR);
    }
    nts = dimension/n_pulse;
    memp.create(nts, &help, &kurtosis_help);
}

/**
 * Function to find the total number of orbitals corresponding to a problem
 * @param readin the problem to be solved
 * @return the number of different orbitals
 */
int MxElem::CalcTotalNumberOfOrbitals(input & readin){

   int counter = 0;
   int Ncores = readin.gNcores();

   for (int cnt=0; cnt<Ncores; cnt++){

      Gauss * atom = readin.gGaussInfo(cnt);
      int Ntypes = atom->gNtypes();
      
      for (int cnt2=0; cnt2<Ntypes; cnt2++){

         char type = atom->gtype(cnt2);
         int L = GetLofType(type);
         counter += ((L+1)*(L+2))/2;

      }

   }

   return counter;

}

int input_stack::peek_char()
{
  while (current_input != 0) {
    int c = current_input->peek();
    if (c != EOF)
      return c;
    if (current_input->next == 0)
      return EOF;
    input *tem = current_input;
    current_input = current_input->next;
    delete tem;
  }
  return EOF;
}

void rand_planck::init( input & args )
{
    int n_x( 0 );
    args.find_key( "planck_res", n_x, 100 );

    x_vec.clear(  );    
    const double dx = 1. / ( n_x - 1 );
    const int max_x = 20;	// "Magic" -> precision lim
    for( int i = 0; i < max_x * n_x; ++ i )
	x_vec.push_back( max_x - i * dx );
    t_vec.push_back( 0 );	// Just let it go.

    this->integrate(  );
    return;
}

void copy(output& dest, input& src)
{
    uint8_t buffer[BUFFER_SIZE];
    stream::len total_written = 0;
    stream::len r;
    do {
        r = src.try_read(buffer, sizeof(buffer));
        if (r == 0) break;
        stream::len w = dest.try_write(buffer, r);
        total_written += w;
        if (w < r) {
            // Did not write the full buffer
            throw incomplete_write(total_written);
        }
    } while (r == sizeof(buffer));
    return;
}

/**
 * Copy constructor for the input class
 * @param lisa the input to be copied
 */
input::input(input & lisa){

   initelements();
   Charge = lisa.gCharge();
   RotationSymm = lisa.gRotationSymm();
   basisset = lisa.gbasisset();
   Ncores = lisa.gNcores();
   cores = new int[Ncores];
   vectors = new R*[Ncores];
   GaussInfo = new Gauss*[Ncores];
   for (int cnt=0; cnt<Ncores; cnt++){
      cores[cnt] = lisa.gcore(cnt);
      vectors[cnt] = new R(*(lisa.gvector(cnt)));
      GaussInfo[cnt] = new Gauss(*(lisa.gGaussInfo(cnt)));
   }

}

void jones_optical::postprocess(input& invals){
#ifdef gen_t_dat
    func_dat=fopen("python/grad_data2.out", "w");
    func_score = fopen("python/score_data2.out", "w");
#else
    func_dat=fopen("python/grad_data.out", "w");
    func_score = fopen("python/score_data.out", "w");
#endif
    stable_spectral_pde_1d_tmpl<comp>::postprocess(invals);

    if(dimension%2){
        err("System jones_optical requires an even dimension", "jones_optical::postprocess",
                "system/jones_optical.cpp", FATAL_ERROR);
    }
    int num_segments;
    invals.retrieve(num_segments, "num_jones_segments", this);
    if(num_segments < 0){
        err("Number of jones segments must be greater than or equal to zero",
                "jones_optical::postprocess", "system/jones_optical.cpp", FATAL_ERROR);
    }
    invals.retrieve(jones_int_dist, "jones_int_dist", this);
    if(jones_int_dist<0){
        err("The distance between jones segments, jones_int_dist,  must be greater than or equal to zero",
                "jones_optical::postprocess", "system/jones_optical.cpp", FATAL_ERROR);
    }

    memp.add(dimension, &nvec1);
    memp.add(nts, &help, &t, &kurtosis_help, &phold, &nvec2);
    double dt = 60.0/nts;
    gaussian_noise(ucur, dimension, 0, 0.2); 
    for(size_t i = 0; i < nts; i++){
        t[i] = dt*(i-nts/2.0);
    }
    for(size_t i = 0; i < nts; i++){
        ucur[i] = ucur[i+nts] = 1.00/cosh(t[i]/2.0);
        help[i] = _real(ucur[i]);
        nvec1[i] = ucur[i];
    }
    for(size_t i = 0; i < num_pulses; i++){
        fft(nvec1 + i*nts, nvec1 +1*nts, nts);
    }
    //generate variables for the jones matrices
    //create jones matrices
    std::string name_base = "jones_system_vars";
    std::string mat_base = "jones_system_matrices";

    for(int i = 0; i < num_segments; i++){
        std::vector<std::shared_ptr<variable> > vv(4);
        for(auto& val : vv){
            val = std::make_shared<variable>();
            val->setname(get_unique_name(name_base));
            val->holder=holder;
            val->parse("0.1");
#ifdef gen_t_dat
            val->set(0*2*3.1415*(rand()*1.0/RAND_MAX));
#else
            val->set(0);
#endif
            invals.insert_item(val);
            cont->addvar(val);
        }
        std::shared_ptr<jones_matrix> m = std::make_shared<jones_matrix>(get_unique_name(mat_base), i, holder);
        invals.insert_item(m);
        m->setup(vv);
        jones_matrices.push_back(m);
    }
}

/**
 * Initialise the matrix elements
 * @param readin the problem to be solved
 */
void MxElem::Init(input & readin){

   MxElemFiller filler(readin);
   int count1, count2, count3, count4, start, start2, start3;
   double OneBodyElement;

   count1 = -1;

   for (int i=0; i<readin.gNcores(); i++){
      Gauss * first = readin.gGaussInfo(i);

      for (int k=0; k<(first->gNtypes()); k++){
         int L1 = GetLofType(first->gtype(k));

         // loop over different [n1x,n1y,n1z] contributions with n1x+n1y+n1z=L1;
         for (int n1x=L1; n1x>=0; n1x--){
            for (int n1y=L1-n1x; n1y>=0; n1y--){
               int n1z = L1-n1x-n1y;
               count1++;

               count2 = -1;

               for (int j=i; j<readin.gNcores(); j++){
                  Gauss * second = readin.gGaussInfo(j);
                  start = 0;
                  if (i==j) start=k;

                  for (int l=start; l<(second->gNtypes()); l++){
                     int L2 = GetLofType(second->gtype(l));

                     // loop over different [n2x,n2y,n2z] contributions with n2x+n2y+n2z=L2;
                     for (int n2x=L2; n2x>=0; n2x--){
                        for (int n2y=L2-n2x; n2y>=0; n2y--){
                           int n2z = L2-n2x-n2y;

                           if ((i==j)&&(l==k)){
                              if (n2x>n1x){
                                 n2x=n1x;
                                 n2y=n1y;
                                 n2z=n1z;
                              } else {
                                 if (n2x==n1x){
                                    if (n2y>n1y){
                                       n2y=n1y;
                                       n2z=n1z;
                                    }
                                 }
                              }
                           }

                           count2++;
                           //If they're centered on the same atom and n1i+n2i odd for i=x,y or z -> Overlap & KE 0.0
                           if ((i==j)&&((((n1x+n2x)%2)!=0)||(((n1y+n2y)%2)!=0)||(((n1z+n2z)%2)!=0))){
                              setSoverlap(count1,count1+count2,0.0);
                              setKEseparate(count1,count1+count2,0.0);
                              OneBodyElement = 0.0;
                           } else {
                              setSoverlap(count1,count1+count2,filler.Overlap(i,k,n1x,n1y,n1z,j,l,n2x,n2y,n2z));
                              OneBodyElement = filler.KE(i,k,n1x,n1y,n1z,j,l,n2x,n2y,n2z);
                              setKEseparate(count1,count1+count2,OneBodyElement);
                           }
                           
                           for (int Ncore=0; Ncore<readin.gNcores(); Ncore++)
                              OneBodyElement += filler.ElNucl(i,k,n1x,n1y,n1z,j,l,n2x,n2y,n2z,Ncore);

                           setTelem(count1,count1+count2,OneBodyElement);

                           count3 = -1;

                           for (int m=i; m<readin.gNcores(); m++){
                              Gauss * third = readin.gGaussInfo(m);
                              start2 = 0;
                              if (i==m) start2=k;

                              for (int n=start2; n<(third->gNtypes()); n++){
                                 int L3 = GetLofType(third->gtype(n));
                                 for (int n3x=L3; n3x>=0; n3x--){
                                    for (int n3y=L3-n3x; n3y>=0; n3y--){
                                       int n3z = L3-n3x-n3y;

                                       if ((i==m)&&(n==k)){
                                          if (n3x>n1x){
                                             n3x=n1x;
                                             n3y=n1y;
                                             n3z=n1z;
                                          } else {
                                             if (n3x==n1x){
                                                if (n3y>n1y){
                                                   n3y=n1y;
                                                   n3z=n1z;
                                                }
                                             }
                                          }
                                       }

                                       count3++;
//.........这里部分代码省略.........

/*

In the diagram below Ui represents the interior state. This should be represnted by (*this). Ug1 and Ug2 are the first and second layers of ghost cells respectively. They
are the output of the function. The input "layer" determines whether Ug1 or Ug2 is output for BCs that use extrapolation. Reflected boundaries (slipWall, viscousWall) do 
not use the layer input. Instead Ui+1 must be specified as (*this) to output Ug2.

____________________|___________
|         |         |          |
|         |         |          |
|   Ug2   |   Ug1   |   Ui     |
|         |         |          |
|         |         |          |
|_________|_________|__________|
                    |
             boundary face

Currently the following boundary conditions are supported: slipWall, viscousWall, characteristic, stagnationInlet, pressureOutlet, subsonicInflow, subsonicOutflow,
supersonicInflow, supersonicOutflow

*/
primVars primVars::GetGhostState( const string &bcType, const vector3d<double> &areaVec, const string &surf, const input &inputVars, const idealGas &eqnState, const int layer )const{
  //bcType -- type of boundary condition to supply ghost cell for
  //areaVec -- area vector of boundary face
  //surf -- i, j, k surface of boundary
  //inputVar -- all input variables
  //eqnState -- equation of state
  //layer -- layer of ghost cell to return (first (closest) or second (farthest))
  //the instance of primVars being acted upon should be the interior cell bordering the boundary

  primVars ghostState = (*this);  //set ghost state equal to boundary state to start

  //check to see that ghost layer corresponds to allowable number
  if ( !(layer == 1 || layer == 2) ){
    cerr << "ERROR: Error in primVars::GetGhostState. Requesting ghost state at a ghost layer " << layer << ". Please choose either 1 or 2" << endl;
    exit(0);
  }

  //normalize area vector (should always point out of domain)
  vector3d<double> normArea;
  if (surf == "il" || surf == "jl" || surf == "kl"){
    normArea = -1.0 * areaVec / areaVec.Mag(); //at lower surface normal should point out of domain for ghost cell calculation
  }
  else if (surf == "iu" || surf == "ju" || surf == "ku"){
    normArea = areaVec / areaVec.Mag(); 
  }

  double normVelCellCenter = 0;

  //slip wall boundary condition ----------------------------------------------------------------------------------------------------------------
  //slipWall is implemented as a reflection of the interior states so that there is no mass flow across the boundary.
  if (bcType == "slipWall"){             //for slip wall state should be reflected across boundary face, density and pressure stay equal to the boundary cell
    vector3d<double> stateVel = (*this).Velocity();
    normVelCellCenter = stateVel.DotProd(normArea);

    //for a slip wall the velocity of the boundary cell center is reflected across the boundary face to get the velocity at the ghost cell center
    vector3d<double> ghostVel (stateVel.X() - 2.0 * normArea.X() * normVelCellCenter,
			       stateVel.Y() - 2.0 * normArea.Y() * normVelCellCenter,
			       stateVel.Z() - 2.0 * normArea.Z() * normVelCellCenter);

    ghostState.data[1] = ghostVel.X();
    ghostState.data[2] = ghostVel.Y();
    ghostState.data[3] = ghostVel.Z();

    //numerical BCs for rho and pressure, same as boundary state

  }
  //viscous wall boundary condition ----------------------------------------------------------------------------------------------------------------
  //viscous wall uses the interior density and pressure, but flips the sign on the velocity so that the velocity at the boundary is 0.
  else if (bcType == "viscousWall"){             //for viscous wall velocity at face should be 0.0, density and pressure stay equal to the boundary cell
    vector3d<double> stateVel = (*this).Velocity();

    //ghost cell velocity at cell center is set to opposite of velocity at boundary cell center so that velocity at face will be zero
    vector3d<double> ghostVel (-1.0 * stateVel.X(),
			       -1.0 * stateVel.Y(),
			       -1.0 * stateVel.Z() );

    ghostState.data[1] = ghostVel.X();
    ghostState.data[2] = ghostVel.Y();
    ghostState.data[3] = ghostVel.Z();

    //numerical BCs for rho and pressure, same as boundary state

  }
  //subsonic inflow boundary condition --------------------------------------------------------------------------------------------------------------
  //this boundary condition enforces density and velocity as freestream inputs (constant) and extrapolates from the interior state to get pressure
  //this is a primative implementation, stagnationInlet or characteristic are better options
  else if (bcType == "subsonicInflow"){     //set velocity and density to freestream values
    double sos = eqnState.GetSoS( inputVars.PRef(), inputVars.RRef() );
    vector3d<double> ghostVel = inputVars.VelRef() / sos;   //nondimensionalize velocity

    ghostState.data[0] = 1.0;
    ghostState.data[1] = ghostVel.X();
    ghostState.data[2] = ghostVel.Y();
    ghostState.data[3] = ghostVel.Z();

    //numerical bc for pressure, same as boundary state

    if (layer == 2){ //extrapolate to get ghost state at 2nd layer
      ghostState = 2.0 * ghostState - (*this);
    }
//.........这里部分代码省略.........

/*!
 * This function initializes the internal memory and variables in the initiator class
 * When called, all parameters returned by dependencies are promised to have been processed
 * at the appropriate level of detail
 * \sa item_dim::postprocess
 */
void integrator::postprocess(input& dat){
    item_dim::postprocess(dat); 
    dat.retrieve(rh_val, "rhs", this);
}