void InterpolateBicubic::set_table( const std::vector<Value>& ff ){
  plumed_assert( getNumberOfSplinePoints()==ff.size() ); 
  plumed_assert( ff[0].getNumberOfDerivatives()==2 );

  dcross=0.0; unsigned iplus, iminus;
  for(unsigned i=1;i<np[0]-1;++i){
      iplus=(i+1)*stride[0]; iminus=(i-1)*stride[0];
      for(unsigned j=1;j<np[1]-1;++j){
          dcross(i,j) = ( ff[iplus+j+1].get() + ff[iminus+j-1].get() - ff[iplus+j-1].get() - ff[iminus+j+1].get() ) /
                          getCrossTermDenominator( i, j );
      }
  } 

  double d1, d2; Matrix<double> tc(4,4);
  std::vector<double> y(4), dy1(4), dy2(4), d2y12(4);

  unsigned pij=0; unsigned ipos;
  for (unsigned i=0;i<np[0]-1;++i){
      ipos=i*stride[0]; d1 = getPointSpacing( 0, i );    
      for (unsigned j=0; j<np[1]-1;++j){
         d2 = getPointSpacing( 1, j );                   
         y[0] = ff[ipos+j].get(); y[1] = ff[ipos+stride[0]+j].get(); y[2] = ff[ipos+stride[0]+j+1].get(); y[3] = ff[ipos+j+1].get();
         dy1[0] = ff[ipos+j].getDerivative(0); dy1[1] = ff[ipos+stride[0]+j].getDerivative(0); 
         dy1[2] = ff[ipos+stride[0]+j+1].getDerivative(0); dy1[3] = ff[ipos+j+1].getDerivative(0);
         dy2[0] = ff[ipos+j].getDerivative(1); dy2[1] = ff[ipos+stride[0]+j].getDerivative(1); 
         dy2[2] = ff[ipos+stride[0]+j+1].getDerivative(1); dy2[3] = ff[ipos+j+1].getDerivative(1);
         d2y12[0] = dcross( i, j ); d2y12[1] = dcross( i+1, j ); d2y12[2] = dcross( i+1, j+1 ); d2y12[3] = dcross( i, j+1 );
         IBicCoeff( y, dy1, dy2, d2y12, d1, d2, tc);

         pij=( ipos+j )*16;
         for(unsigned k=0; k<4; ++k){ for(unsigned n=0; n<4; ++n){ clist[pij++]=tc(k,n); } }
      }
  }
}

void ActionWithInputVessel::readArgument( const std::string& type ) {
  std::string mlab;
  if( keywords.exists("DATA") && type!="grid" ) parse("DATA",mlab);
  ActionWithVessel* mves= plumed.getActionSet().selectWithLabel<ActionWithVessel*>(mlab);
  if(!mves) error("action labelled " +  mlab + " does not exist or does not have vessels");
  addDependency(mves);

  ActionWithValue* aval=dynamic_cast<ActionWithValue*>( this );
  if(aval) {
    if( aval->checkNumericalDerivatives() ) {
      ActionWithValue* aval2=dynamic_cast<ActionWithValue*>( mves );
      plumed_assert( aval2 ); aval2->useNumericalDerivatives();
    }
  }

  if( type=="bridge" ) {
    ActionWithVessel* aves=dynamic_cast<ActionWithVessel*>( this );
    plumed_assert(aves); myBridgeVessel = mves->addBridgingVessel( aves );
    arguments = dynamic_cast<Vessel*>( myBridgeVessel );
  } else  if( type=="store" ) {
    arguments = dynamic_cast<Vessel*>( mves->buildDataStashes( NULL ) );
  } else {
    plumed_error();
  }
}

void ActionWithVessel::addTaskToList( const unsigned& taskCode ){
  plumed_assert( functions.size()==0 );  // We cannot add more tasks after vessels added
  indexOfTaskInFullList.push_back( fullTaskList.size() );
  fullTaskList.push_back( taskCode ); partialTaskList.push_back( taskCode ); 
  taskFlags.push_back(0); nactive_tasks = fullTaskList.size();
  plumed_assert( partialTaskList.size()==nactive_tasks && indexOfTaskInFullList.size()==nactive_tasks && taskFlags.size()==nactive_tasks );
}

void PlumedMain::readInputWords(const std::vector<std::string> & words) {
  plumed_assert(initialized);
  if(words.empty())return;
  else if(words[0]=="ENDPLUMED") return;
  else if(words[0]=="_SET_SUFFIX") {
    plumed_assert(words.size()==2);
    setSuffix(words[1]);
  } else {
    std::vector<std::string> interpreted(words);
    Tools::interpretLabel(interpreted);
    Action* action=actionRegister().create(ActionOptions(*this,interpreted));
    if(!action) {
      log<<"ERROR\n";
      log<<"I cannot understand line:";
      for(unsigned i=0; i<interpreted.size(); ++i) log<<" "<<interpreted[i];
      log<<"\n";
      log.flush();
      plumed_merror("I cannot understand line " + interpreted[0] + " " + interpreted[1]);
    };
    action->checkRead();
    actionSet.push_back(action);
  };

  pilots=actionSet.select<ActionPilot*>();
}

void LinkCells::buildCellLists( const std::vector<Vector>& pos, const std::vector<unsigned>& indices, const Pbc& pbc ){
  plumed_assert( cutoffwasset && pos.size()==indices.size() );

  // Must be able to check that pbcs are not nonsensical in some way?? -- GAT

  // Setup the pbc object by copying it from action
  mypbc.setBox( pbc.getBox() );

  // Setup the lists
  if( pos.size()!=allcells.size() ){ 
    allcells.resize( pos.size() ); lcell_lists.resize( pos.size() ); 
  }

   {
// This is the reciprocal lattice
// notice that reciprocal.getRow(0) is a vector that is orthogonal to b and c
// This allows to use linked cells in non orthorhomic boxes
     Tensor reciprocal(transpose(mypbc.getInvBox()));
     ncells[0] = std::floor( 1.0/ reciprocal.getRow(0).modulo() / link_cutoff );
     if( ncells[0]==0 ) ncells[0]=1;
     ncells[1] = std::floor( 1.0/ reciprocal.getRow(1).modulo() / link_cutoff );
     if( ncells[1]==0 ) ncells[1]=1;
     ncells[2] = std::floor( 1.0/ reciprocal.getRow(2).modulo() / link_cutoff );
     if( ncells[2]==0 ) ncells[2]=1;
  }
  // Setup the strides
  nstride[0]=1; nstride[1]=ncells[0]; nstride[2]=ncells[0]*ncells[1];

  // Setup the storage for link cells
  unsigned ncellstot=ncells[0]*ncells[1]*ncells[2];
  if( lcell_tots.size()!=ncellstot ){
      lcell_tots.resize( ncellstot ); lcell_starts.resize( ncellstot );
  }
  // Clear nlcells
  for(unsigned i=0;i<ncellstot;++i) lcell_tots[i]=0;
  // Clear allcells
  allcells.assign( allcells.size(), 0 );

  // Find out what cell everyone is in
  unsigned rank=comm.Get_rank(), size=comm.Get_size();
  for(unsigned i=rank;i<pos.size();i+=size){
      allcells[i]=findCell( pos[i] );
      lcell_tots[allcells[i]]++;
  }
  // And gather all this information on every node
  comm.Sum( allcells ); comm.Sum( lcell_tots );

  // Now prepare the link cell lists
  unsigned tot=0;
  for(unsigned i=0;i<lcell_tots.size();++i){ lcell_starts[i]=tot; tot+=lcell_tots[i]; lcell_tots[i]=0; }
  plumed_assert( tot==pos.size() );

  // And setup the link cells properly
  for(unsigned j=0;j<pos.size();++j){
      unsigned myind = lcell_starts[ allcells[j] ] + lcell_tots[ allcells[j] ];
      lcell_lists[ myind ] = indices[j];
      lcell_tots[allcells[j]]++;
  }
}

void BridgeVessel::completeNumericalDerivatives(){
  unsigned nextra = myOutputAction->getNumberOfDerivatives() - getAction()->getNumberOfDerivatives();
  Matrix<double> tmpder( myOutputValues->getNumberOfComponents(), nextra );
  ActionWithVessel* vval=dynamic_cast<ActionWithVessel*>( myOutputAction );
  for(unsigned i=0;i<nextra;++i){
      vval->bridgeVariable=i; getAction()->calculate();
      for(int j=0;j<myOutputValues->getNumberOfComponents();++j) tmpder(j,i) = myOutputValues->getOutputQuantity(j);
  }
  vval->bridgeVariable=nextra; getAction()->calculate(); 
  plumed_assert( inum==mynumerical_values.size() ); inum=0;  // Reset inum now that we have finished calling calculate
  std::vector<double> base( myOutputValues->getNumberOfComponents() );
  for(int j=0;j<myOutputValues->getNumberOfComponents();++j) base[j] = myOutputValues->getOutputQuantity(j);

  const double delta=sqrt(epsilon);
  ActionAtomistic* aa=dynamic_cast<ActionAtomistic*>( getAction() );
  unsigned nvals=myOutputValues->getNumberOfComponents();
  for(unsigned j=0;j<nvals;++j) ( myOutputValues->copyOutput(j) )->clearDerivatives();   

  if( aa ){
      ActionWithArguments* aarg=dynamic_cast<ActionWithArguments*>( getAction() );
      plumed_assert( !aarg ); Tensor box=aa->getBox(); 
      unsigned natoms=aa->getNumberOfAtoms();
      for(unsigned j=0;j<nvals;++j){
          double ref=( myOutputValues->copyOutput(j) )->get();
          if( ( myOutputValues->copyOutput(j) )->getNumberOfDerivatives()>0 ){
              for(unsigned i=0;i<3*natoms;++i){
                  double d=( mynumerical_values[i*nvals+j] - ref)/delta;
                  ( myOutputValues->copyOutput(j) )->addDerivative(i,d);
              }
              Tensor virial;
              for(int i=0;i<3;i++) for(int k=0;k<3;k++){
                 virial(i,k)=( mynumerical_values[ nvals*(3*natoms + 3*i + k) + j ]-ref)/delta;
              }
              virial=-matmul(box.transpose(),virial);
              for(int i=0;i<3;i++) for(int k=0;k<3;k++) ( myOutputValues->copyOutput(j) )->addDerivative(3*natoms+3*k+i,virial(k,i));
          }
      }
  } else {
      plumed_merror("not implemented or tested yet");
//      unsigned nder=myOutputAction->getNumberOfDerivatives();
//      for(unsigned j=0;j<nvals;++j){
//          double ref=( myOutputValues->copyOutput(j) )->get();
//              for(unsigned i=0;i<nder;++i){
//                  double d=( mynumerical_values[i*nvals+j] - ref)/delta;
//                  ( myOutputValues->copyOutput(j) )->addDerivative(i,d);
//              }
//          }
//      }
  }
  // Add the derivatives wrt to the local quantities we are working with
  for(unsigned j=0;j<nvals;++j){
     unsigned k=0;
     for(unsigned i=getAction()->getNumberOfDerivatives();i<myOutputAction->getNumberOfDerivatives();++i){
        ( myOutputValues->copyOutput(j) )->addDerivative( i, (tmpder(j,k)-base[j])/sqrt(epsilon) ); k++;
     }
  }
}

MultiReferenceBase::MultiReferenceBase( const std::string& type, const bool& checksoff ):
  wasSet(false),
  skipchecks(checksoff),
  mtype(type)
{
  if(checksoff) plumed_assert( mtype.length()==0 );
}

ManyRestraintsBase::ManyRestraintsBase(const ActionOptions& ao):
    Action(ao),
    ActionWithValue(ao),
    ActionPilot(ao),
    ActionWithVessel(ao),
    ActionWithInputVessel(ao)
{
    // Read in the vessel we are action on
    readArgument("bridge");
    aves=dynamic_cast<ActionWithVessel*>( getDependencies()[0] );

    plumed_assert( getDependencies().size()==1 && aves );
    log.printf("  adding restraints on variables calculated by %s action with label %s\n",
               aves->getName().c_str(),aves->getLabel().c_str());

    // Add a task list in order to avoid problems
    for(unsigned i=0; i<aves->getFullNumberOfTasks(); ++i) addTaskToList( aves->getTaskCode(i) );
    // And turn on the derivatives (note problems here because of ActionWithValue)
    turnOnDerivatives();
    needsDerivatives();

    // Now create the vessel
    std::string fake_input="LABEL=bias";
    addVessel( "SUM", fake_input, 0 );
    readVesselKeywords();
}

OFile& OFile::rewind() {
// we use here "hard" rewind, which means close/reopen
// the reason is that normal rewind does not work when in append mode
// moreover, we can take a backup of the file
  plumed_assert(fp);
  clearFields();
  if(gzfp) {
#ifdef __PLUMED_HAS_ZLIB
    gzclose((gzFile)gzfp);
#endif
  } else fclose(fp);
  if(!comm || comm->Get_rank()==0) {
    std::string fname=this->path;
    size_t found=fname.find_last_of("/\\");
    std::string directory=fname.substr(0,found+1);
    std::string file=fname.substr(found+1);
    std::string backup=directory+backstring +".last."+file;
    int check=rename(fname.c_str(),backup.c_str());
    plumed_massert(check==0,"renaming "+fname+" into "+backup+" failed for reason: "+strerror(errno));
  }
  if(gzfp) {
#ifdef __PLUMED_HAS_ZLIB
    gzfp=(void*)gzopen(const_cast<char*>(this->path.c_str()),"w9");
#endif
  } else fp=std::fopen(const_cast<char*>(path.c_str()),"w");
  return *this;
}

void Stopwatch::Watch::pause() {
  plumed_assert(running>0);
  running--;
  if(running!=0) return;
  auto t=std::chrono::duration_cast<std::chrono::nanoseconds>(std::chrono::high_resolution_clock::now()-lastStart);
  lap+=t.count();
}

void ActionWithVessel::readVesselKeywords(){
  // Set maxderivatives if it is too big
  if( maxderivatives>getNumberOfDerivatives() ) maxderivatives=getNumberOfDerivatives();

  // Loop over all keywords find the vessels and create appropriate functions
  for(unsigned i=0;i<keywords.size();++i){
      std::string thiskey,input; thiskey=keywords.getKeyword(i);
      // Check if this is a key for a vessel
      if( vesselRegister().check(thiskey) ){
          plumed_assert( keywords.style(thiskey,"vessel") );
          bool dothis=false; parseFlag(thiskey,dothis);
          if(dothis) addVessel( thiskey, input );

          parse(thiskey,input);
          if(input.size()!=0){ 
                addVessel( thiskey, input );
          } else {
             for(unsigned i=1;;++i){
                if( !parseNumbered(thiskey,i,input) ) break;
                std::string ss; Tools::convert(i,ss);
                addVessel( thiskey, input, i ); 
                input.clear();
             } 
          }
      }
  }

  // Make sure all vessels have had been resized at start
  if( functions.size()>0 ) resizeFunctions();
}

void Action::registerKeywords( Keywords& keys ) {
  plumed_assert( keys.size()==0 );
  keys.add( "hidden", "LABEL", "a label for the action so that its output can be referenced in the input to other actions.  Actions with scalar output are referenced using their label only.  Actions with vector output must have a separate label for every component.  Individual componets are then refered to using label.component" );
  keys.reserve("optional","UPDATE_FROM","Only update this action from this time");
  keys.reserve("optional","UPDATE_UNTIL","Only update this action until this time");
  keys.reserve("optional","RESTART","allows per-action setting of restart (YES/NO/AUTO)");
}

void MultiColvarBase::activateIndexes( const unsigned& istart, const unsigned& number, const std::vector<unsigned>& indexes ){
  plumed_assert( number>0 );
  for(unsigned i=0;i<number-9;i+=3){
      plumed_dbg_assert( indexes[istart+i]%3==0 ); unsigned iatom=indexes[istart+i]/3; 
      atoms_with_derivatives.activate( iatom ); 
  }
}

void ClusterWithSurface::retrieveAtomsInCluster( const unsigned& clust, std::vector<unsigned>& myatoms ) const {
  std::vector<unsigned> tmpat; myclusters->retrieveAtomsInCluster( clust, tmpat );

  // Prevent double counting
  std::vector<bool> incluster( getNumberOfNodes(), false );
  for(unsigned i=0;i<tmpat.size();++i) incluster[tmpat[i]]=true;

  // Find the atoms in the the clusters
  std::vector<bool> surface_atom( getNumberOfNodes(), false ); 
  for(unsigned i=0;i<tmpat.size();++i){
      for(unsigned j=0;j<getNumberOfNodes();++j){
         if( incluster[j] ) continue;
         double dist2=getSeparation( getPosition(tmpat[i]), getPosition(j) ).modulo2();
         if( dist2<rcut_surf2 ){ surface_atom[j]=true; }
      }
  }
  unsigned nsurf_at=0; 
  for(unsigned j=0;j<getNumberOfNodes();++j){
     if( surface_atom[j] ) nsurf_at++; 
  }
  myatoms.resize( nsurf_at + tmpat.size() );
  for(unsigned i=0;i<tmpat.size();++i) myatoms[i]=tmpat[i];
  unsigned nn=tmpat.size();
  for(unsigned j=0;j<getNumberOfNodes();++j){
      if( surface_atom[j] ){ myatoms[nn]=j; nn++; }
  }
  plumed_assert( nn==myatoms.size() );
}

void HistogramBead::generateBins( const std::string& params, std::vector<std::string>& bins ) {
  std::vector<std::string> data=Tools::getWords(params);
  plumed_massert(data.size()>=1,"There is no input for this keyword");

  std::string name=data[0];

  unsigned nbins; std::vector<double> range(2); std::string smear;
  bool found_nb=Tools::parse(data,"NBINS",nbins);
  plumed_massert(found_nb,"Number of bins in histogram not found");
  bool found_r=Tools::parse(data,"LOWER",range[0]);
  plumed_massert(found_r,"Lower bound for histogram not specified");
  found_r=Tools::parse(data,"UPPER",range[1]);
  plumed_massert(found_r,"Upper bound for histogram not specified");
  plumed_massert(range[0]<range[1],"Range specification is dubious");
  bool found_b=Tools::parse(data,"SMEAR",smear);
  if(!found_b) { Tools::convert(0.5,smear); }

  std::string lb,ub; double delr = ( range[1]-range[0] ) / static_cast<double>( nbins );
  for(unsigned i=0; i<nbins; ++i) {
    Tools::convert( range[0]+i*delr, lb );
    Tools::convert( range[0]+(i+1)*delr, ub );
    bins.push_back( name + " " +  "LOWER=" + lb + " " + "UPPER=" + ub + " " + "SMEAR=" + smear );
  }
  plumed_assert(bins.size()==nbins);
}

void ReferenceArguments::setReferenceArguments( const std::vector<double>& arg_vals, const std::vector<double>& sigma ) {
  moveReferenceArguments( arg_vals );

  if( hasmetric ) {
    unsigned k=0;
    for(unsigned i=0; i<reference_args.size(); ++i) {
      for(unsigned j=i; j<reference_args.size(); ++j) {
        metric(i,j)=metric(j,i)=sigma[k]; k++;
      }
    }
    plumed_assert( k==sigma.size() );
  } else {
    plumed_assert( reference_args.size()==sigma.size() );
    for(unsigned i=0; i<reference_args.size(); ++i) weights[i]=sigma[i];
  }
}

void MultiReferenceBase::getAtomAndArgumentRequirements( std::vector<AtomNumber>& atoms, std::vector<std::string>& args ) {
  plumed_assert( atoms.size()==0 && args.size()==0 );
  for(unsigned i=0; i<frames.size(); ++i) {
    frames[i]->getAtomRequests( atoms );
    frames[i]->getArgumentRequests( args );
  }
}

StoreCentralAtomsVessel::StoreCentralAtomsVessel( const vesselbase::VesselOptions& da ):
StoreDataVessel(da),
tmpdf(3)
{
  mycolv=dynamic_cast<MultiColvarBase*>( getAction() );
  plumed_assert( mycolv ); completeSetup( mycolv->getCentralAtomElementIndex(), 3 );
}

void TargetDist::read( const std::vector<double>& targ, std::vector<Value*> ar ) {
  plumed_assert( targ.size()==ar.size() );

  target.resize( ar.size() ); args.resize( ar.size() );
  log.printf("  distance from this point in cv space : ");
  for(unsigned i=0; i<target.size(); ++i) { log.printf("%f ", targ[i]); target[i]=targ[i]; args[i]=ar[i]; }
  log.printf("\n");
}

void MDAtomsTyped<T>::setp(void*pp){
  T*p=static_cast<T*>(pp);
  plumed_assert(stride==0 || stride==3);
  px=p;
  py=p+1;
  pz=p+2;
  stride=3;
}