void assign(state_type state, tree_type& tree)
    {
      features_ += *state.feature_vector();

      switch (state.operation().operation()) {
      case operation_type::AXIOM:
	break;
      case operation_type::FINAL:
      case operation_type::IDLE:
	assign(state.derivation(), tree);
	break;
      case operation_type::UNARY:
	tree.label_ = state.label();
	tree.antecedent_.resize(1);
	assign(state.derivation(), tree.antecedent_.front());
	break;
      case operation_type::SHIFT:
	tree.label_ = state.label();
	tree.antecedent_ = tree_type::antecedent_type(1, state.head());
	break;
      case operation_type::REDUCE:
      case operation_type::REDUCE_LEFT:
      case operation_type::REDUCE_RIGHT:
	tree.label_ = state.label();
	tree.antecedent_.resize(2);
	assign(state.reduced(), tree.antecedent_.front());
	assign(state.derivation(), tree.antecedent_.back());
	break;
      }
    }

int saveHdf5(state_type &OUT, state_type &TIME){
    hid_t   hdf_file,hdf_group,hdf_data,dataspace_id;
    herr_t  status;

    fprintf(stdout,"Writing file %s ...",allparams.outfilename);
    hdf_file = H5Fcreate(allparams.outfilename,H5F_ACC_TRUNC,H5P_DEFAULT,H5P_DEFAULT);
    if (hdf_file < 0){
        return -1;
    }

    /*
    if ( (hdf_group=H5Gopen2(hdf_file,"/",H5P_DEFAULT)) < 0){
        H5Fclose(hdf_file);
        return -1;
    }*/
    
    /* Write particle positions and velocities.
     * Ordered in chunk where first Nstep+1 lines correspond to particle 1,
     * step Nstep+1 chunk correspond to particle 2 etc. */
    std::cout << "Writing positions and velocities\n";
    hsize_t dims[1]={OUT.size()};
    dataspace_id=H5Screate_simple(1,dims,NULL);
    if ( (hdf_data=H5Dcreate2(hdf_file,"x",H5T_NATIVE_DOUBLE,dataspace_id,H5P_DEFAULT,H5P_DEFAULT,H5P_DEFAULT)) < 0){
        H5Dclose(hdf_data);
        return -1;
    }
    status=H5Dwrite(hdf_data, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, &OUT[0]);

    /*Write times*/
    std::cout << "Writing times\n";
    hsize_t dims2[1]={TIME.size()};
    dataspace_id=H5Screate_simple(1,dims2,NULL);
    if ( (hdf_data=H5Dcreate2(hdf_file,"t",H5T_NATIVE_DOUBLE,dataspace_id,H5P_DEFAULT,H5P_DEFAULT,H5P_DEFAULT)) < 0){
        H5Dclose(hdf_data);
        return -1;
    }
    status=H5Dwrite(hdf_data, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, &TIME[0]);

    /*Write no. of components*/
    dims2[0]={1};
    dataspace_id=H5Screate_simple(1,dims2,NULL);
    if ( (hdf_data=H5Dcreate2(hdf_file,"NumComponents",H5T_NATIVE_INT,dataspace_id,H5P_DEFAULT,H5P_DEFAULT,H5P_DEFAULT)) < 0){
        H5Dclose(hdf_data);
        return -1;
    }
    status=H5Dwrite(hdf_data, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, &allparams.NumComponents);


    H5Fclose(hdf_file);
    H5Dclose(hdf_data);
    H5Sclose(dataspace_id);
    fprintf(stdout," file written successfully!\n");
    return 0;
}

double compute_profile(const vector<double>& rate, state_type& profile)
{
  polysome tasep_tmp(&rate);
  if (profile.size() != rate.size()-1) {
    cout<<profile.size()<<" "<<rate.size()<<endl;
    cout<<"profile length not equal to rate vector length-1!"<<endl;
    return -1;
  }
  tasep_tmp.run();
  profile = tasep_tmp.get_Aprob();
  return tasep_tmp.compute_translation_rate();
}

void
FMUIntegrator::integrate( fmiReal step_size, fmiReal dt )
{
	// This vector holds (temporarily) the values of the FMU's continuous states.
	static state_type states( fmu_->nStates() );

	// Get current continuous states.
	fmu_->getContinuousStates( &states.front() );

	// Invoke integration method.
  	stepper_->invokeMethod( this, states, fmu_->getTime(), step_size, dt );
}

 static void resize( state_type &x1 ,
                     const state_type &x2 )
 {
     // allocate required memory
     x1.resize( x2.size() );
     for( size_t i=0 ; i < x2.size() ; ++i )
     {
         x1[i] = dataflow< hpx_resize_2d_action >( find_here() , 
                                                   std::allocate_shared<dvecvec>( std::allocator<dvecvec>() ) , 
                                                   x2[i].get_future().get() );
     }
 }

        void operator()(const state_type& x, state_type& dxdt, const double& t)
        {
            for (state_type::iterator i(dxdt.begin()); i != dxdt.end(); ++i)
            {
                *i = 0.0;
            }
            // XXX
            for (reaction_container_type::const_iterator
                i(reactions_.begin()); i != reactions_.end(); i++)
            {
                // Prepare  state_array of reactants and products that contain amounts of each reactants.
                Ratelaw::state_container_type reactants_states(i->reactants.size());
                Ratelaw::state_container_type products_states(i->products.size());
                Ratelaw::state_container_type::size_type cnt(0);

                for (index_container_type::const_iterator
                    j((*i).reactants.begin()); j != (*i).reactants.end(); ++j, cnt++)
                {
                    reactants_states[cnt] = x[*j];
                }
                cnt = 0;
                for (index_container_type::const_iterator
                    j((*i).products.begin()); j != (*i).products.end(); ++j, cnt++)
                {
                    products_states[cnt] = x[*j];
                }

                double flux;
                // Get pointer of Ratelaw object and call it.
                if (i->ratelaw.expired() || i->ratelaw.lock()->is_available() == false)
                {
                    boost::scoped_ptr<Ratelaw> temporary_ratelaw_obj(new RatelawMassAction(i->k));
                    flux = temporary_ratelaw_obj->deriv_func(reactants_states, products_states, volume_);
                }
                else
                {
                    boost::shared_ptr<Ratelaw> ratelaw = (*i).ratelaw.lock();
                    flux = (*ratelaw).deriv_func(reactants_states, products_states, volume_);
                }

                // Merge each reaction's flux into whole dxdt
                for (index_container_type::const_iterator
                    j((*i).reactants.begin()); j != (*i).reactants.end(); ++j)
                {
                    dxdt[*j] -= flux;
                }
                for (index_container_type::const_iterator
                    j((*i).products.begin()); j != (*i).products.end(); ++j)
                {
                    dxdt[*j] += flux;
                }
            }
        }

void lorenz_with_lyap( const state_type &x , state_type &dxdt , double t )
{
    lorenz()( x , dxdt , t );

    for( size_t l=0 ; l<num_of_lyap ; ++l )
    {
        const double *pert = x.begin() + 3 + l * 3;
        double *dpert = dxdt.begin() + 3 + l * 3;
        dpert[0] = - sigma * pert[0] + 10.0 * pert[1];
        dpert[1] = ( R - x[2] ) * pert[0] - pert[1] - x[0] * pert[2];
        dpert[2] = x[1] * pert[0] + x[0] * pert[1] - b * pert[2];
    }
}

void lattice( const state_type &x , state_type &dxdt , const double /* t */ )
{
    state_type::const_iterator x_begin = x.begin();
    state_type::const_iterator x_end = x.end();
    state_type::iterator dxdt_begin = dxdt.begin();

    x_end--; // stop one before last
    while( x_begin != x_end )
    {
        *(dxdt_begin++) = std::sin( *(x_begin) - *(x_begin++) );
    }
    *dxdt_begin = sin( *x_begin - *(x.begin()) ); // periodic boundary
}

void
FMUIntegrator::operator()( const state_type& x, state_type& dx, fmiReal time )
{
	// if there has been an event, then the integrator shall do nothing
	if ( ! fmu_->getStateEventFlag() ) {
		// Update to current time.
		fmu_->setTime( time );

		// Update to current states.
		fmu_->setContinuousStates( &x.front() );

		// Evaluate derivatives and store them to vector dx.
		fmu_->getDerivatives( &dx.front() );
	}
}

  void operator() (const state_type x, state_type& dxdt, double t) {

    // for input voltage, model a step function at time zero
    double Va = 0;
    if (t > 0.0) {
      Va = vin_;
    }
    Matrix<double, 1, 1> u; u << Va;

    // All other node voltages are determined by odeint through our equations:
    Map<const Matrix<double, 2, 1> > xvec(x.data());
    Map<Matrix<double, 2, 1> > result(dxdt.data());

    result = coeff_ * xvec + input_ * u;
  }

void system_2d( const state_type &q , state_type &dpdt )
{
    // works on shared data, but coupling data is provided as copy
    const size_t N = q.size();
    //state_type dpdt_(N);
    // first row
    dpdt[0] = dataflow< system_first_block_action >( find_here() , q[0] , 
                                                     dataflow< first_row_action >( find_here() , q[1] ) , 
                                                     dpdt[0] , 0 );
    // middle rows
    for( size_t i=1 ; i<N-1 ; i++ )
    {
        dpdt[i] = dataflow< system_center_block_action >( find_here() , q[i] , 
                                                          dataflow< last_row_action >( find_here() , q[i-1] ) , 
                                                          dataflow< first_row_action >( find_here() , q[i+1] ) , 
                                                          dpdt[i] , i );
    }
    dpdt[N-1] = dataflow< system_last_block_action >( find_here() , q[N-1] , 
                                                      dataflow< last_row_action >( find_here() , q[N-2] ) , 
                                                      dpdt[N-1] , N-1);

    // coupling synchronization step
    // dpdt[0] = dataflow< sys_sync1_action >( fing_here() , dpdt_[0] , dpdt_[1] );
    // for( size_t i=1 ; i<N-1 ; i++ )
    // {
    //     dpdt[i] = dataflow< sys_sync2_action >( find_here() , dpdt_[i] , dpdt_[i] , q[i+1] , dpdt[i] );
    // }
    // dpdt_[N-1] = dataflow< system_last_block_action >( find_here() , q[N-1] , q[N-2] , dpdt[N-1] );

}

 std::codecvt_base::result 
 do_in( state_type& state, const char* first1, const char* last1, 
        const char*& next1, wchar_t* first2, wchar_t* last2, 
        wchar_t*& next2 ) const
 {
     using namespace std;
     if (state < 0 || state > 3)
         return codecvt_base::error;
     next1 = first1;
     next2 = first2;
     while (next2 != last2 && next1 != last1) {
         while (next1 != last1) {
             if (state == 0) {
                 if (*next1 < 1 || *next1 > 3)
                     return codecvt_base::error;
                 state = *next1++;
             } else if (state == 1) {
                 *next2++ = (unsigned char) *next1++;
                 state = 0;
                 break;
             } else {
                 if (*next1++ != 0)
                     return codecvt_base::error;
                 --state.val();
             }
         }
     }
     return next2 == last2 ? 
         codecvt_base::ok : 
         codecvt_base::partial;
 }

 std::codecvt_base::result 
 do_out( state_type& state, const wchar_t* first1, const wchar_t* last1,
         const wchar_t*& next1, char* first2, char* last2, 
         char*& next2 ) const
 {
     using namespace std;
     if (state < 0 || state > 3)
         return codecvt_base::error;
     next1 = first1;
     next2 = first2;
     while (next1 != last1 && next2 != last2) {
         while (next2 != last2) {
             if (state == 0) {
                 if (*next1 > integer_traits<unsigned char>::const_max)
                     return codecvt_base::noconv;
                 state = *next1 % 3 + 1;
                 *next2++ = static_cast<char>(state);
             } else if (state == 1) {
                 state = 0;
                 *next2++ = static_cast<unsigned char>(*next1++);
                 break;
             } else {
                 --state.val();
                 *next2++ = 0;
             }
         }
     }
     return next1 == last1 ? 
         codecvt_base::ok : 
         codecvt_base::partial;
 }

typename StackedGatedModel<Z>::State StackedGatedModel<Z>::activate(
        state_type& previous_state,
        const Indexing::Index indices) const {
    State out;
    auto input_vector = this->embedding[indices];

    out.memory       = gate.activate(
        {
            input_vector,
            previous_state.back().hidden
        }
    ).sigmoid();

    input_vector      = input_vector.eltmul_broadcast_colwise(out.memory);

    out.lstm_state = this->stacked_lstm.activate(
        previous_state,
        input_vector);

    out.prediction = MatOps<Z>::softmax_rowwise(
            this->decode(
                input_vector,
                out.lstm_state
            )
        );

    return out;
}

 void operator() (const state_type psit, double t) const
 {
     // If the above line doesn't work for you, sub in the one below.
     m_out << t;        
     for (size_t i=0; i < psit.size(); ++i) m_out << "\t" << psit(i).real() << "\t" << psit(i).imag();
     m_out << "\n";
 }

 void operator()( const state_type &x , state_type &dxdt , const value_type t ) const
 {
     std::cout << "sys start " << dxdt.size() << std::endl;
     dxdt[0] = x[0] + 2 * x[1];
     dxdt[1] = x[1];
     std::cout << "sys done" << std::endl;
 }

    void do_step( System system , state_type& x , time_type t , time_type dt )
    {
        typedef typename odeint::unwrap_reference< System >::type system_type;
        typedef typename odeint::unwrap_reference< typename system_type::first_type >::type deriv_func_type;
        typedef typename odeint::unwrap_reference< typename system_type::second_type >::type jacobi_func_type;
        system_type &sys = system;
        deriv_func_type &deriv_func = sys.first;
        jacobi_func_type &jacobi_func = sys.second;   
        m_resizer.adjust_size( x , detail::bind( &stepper_type::template resize_impl<state_type> , detail::ref( *this ) , detail::_1 ) );
         
        state_type x1 = x;// x1 for storing x(k+1)
        // std::cout << x1[0] << " " << x1[1] << "\n";
        implicit_euler< double > imp_euler;
        imp_euler.do_step(system, x1, t, dt);
        //x(k+1) computed from above line for x(k) = 4*x(k+1) - 3*x(k+2) + 2*dt*x'(k+2)
        
        // std::cout << x1[0] << " " << x1[1] << "\n";
        
        
        for( size_t i=0 ; i<x.size() ; ++i )
            m_pm.m_v[i] = i;
        
        t +=dt;// t = t0 + 2*dt for x(k+2) and x'(k+2)
        deriv_func( x1 , m_dxdt.m_v , t ); //m_dxdt.m_v is x'(k+2)
        m_b.m_v = 2.0 * dt * m_dxdt.m_v;
        jacobi_func( x1 , m_jacobi.m_v  , t );
        m_jacobi.m_v *= dt;
        m_jacobi.m_v -= boost::numeric::ublas::identity_matrix< value_type >( x.size() );
        
        solve( m_b.m_v , m_jacobi.m_v );
        m_x.m_v = (4.0 *x1)- (3.0*x) + m_b.m_v;//x(k) after first Newton iteration
        
        std::cout << "Before loop : " << m_x.m_v[0] << " " << m_x.m_v[1] << " " << m_b.m_v[0] << " " << m_b.m_v[1] << "\n";
        
        
        while( boost::numeric::ublas::norm_2( m_b.m_v ) > m_epsilon )
        {
            deriv_func( m_x.m_v , m_dxdt.m_v , t );
            m_b.m_v = -(4.0 *x1)+ (3.0*x) + m_x.m_v + 2.0 * dt * m_dxdt.m_v;
            solve( m_b.m_v , m_jacobi.m_v );
            m_x.m_v -= m_b.m_v;
            
            std::cout << "Inside loop : " << m_x.m_v[0] << " " << m_x.m_v[1] << " " << m_b.m_v[0] << " " << m_b.m_v[1] << "\n";
        }

        x = m_x.m_v;
    } 

    void do_step( System system , state_type &x , time_type t , time_type dt )
    {
        typedef typename odeint::unwrap_reference< System >::type system_type;
        typedef typename odeint::unwrap_reference< typename system_type::first_type >::type deriv_func_type;
        typedef typename odeint::unwrap_reference< typename system_type::second_type >::type jacobi_func_type;
        system_type &sys = system;
        deriv_func_type &deriv_func = sys.first;
        jacobi_func_type &jacobi_func = sys.second;

        m_resizer.adjust_size( x , detail::bind( &stepper_type::template resize_impl<state_type> , detail::ref( *this ) , detail::_1 ) );

        for( size_t i=0 ; i<x.size() ; ++i )
            m_pm.m_v[i] = i;

        t += dt;

        // apply first Newton step
        deriv_func( x , m_dxdt.m_v , t );

        m_b.m_v = dt * m_dxdt.m_v;

        jacobi_func( x , m_jacobi.m_v  , t );
        m_jacobi.m_v *= dt;
        m_jacobi.m_v -= boost::numeric::ublas::identity_matrix< value_type >( x.size() );

        solve( m_b.m_v , m_jacobi.m_v );

        m_x.m_v = x - m_b.m_v;

        // iterate Newton until some precision is reached
        // ToDo: maybe we should apply only one Newton step -> linear implicit one-step scheme
        while( boost::numeric::ublas::norm_2( m_b.m_v ) > m_epsilon )
        {
            deriv_func( m_x.m_v , m_dxdt.m_v , t );
            m_b.m_v = x - m_x.m_v + dt*m_dxdt.m_v;

            // simplified version, only the first Jacobian is used
            //            jacobi( m_x , m_jacobi , t );
            //            m_jacobi *= dt;
            //            m_jacobi -= boost::numeric::ublas::identity_matrix< value_type >( x.size() );

            solve( m_b.m_v , m_jacobi.m_v );

            m_x.m_v -= m_b.m_v;
        }
        x = m_x.m_v;
    }

 static void resize( state_type &x1 , const state_type &x2 )
 {
     x1.resize( boost::size( x2 ) );
     typename state_type::iterator begin1 = boost::begin( x1 );
     typename state_type::const_iterator begin2 = boost::begin( x2 );
     while( begin1 != boost::end( x1 ) )
         (*begin1++).resize( boost::size( *begin2++ ) );
 }

  // utility function for displaying data.  Applies internal matrices to state
  std::vector<double> state2output(state_type const& x) {
    std::vector<double> result(2);
    Map<const Matrix<double, states, 1> > xvec(x.data());

    Map<Matrix<double, 3, 1> > ovec(result.data());
    ovec = Lred_.transpose() * xvec;  // in theory the input could be involved via Dred_
    return result;
  }