void lcs(seq const & xs, seq const & ys, seq & an_lcs)
{
    members xs_in_lcs(xs.size(), false);
    calculate_lcs(xs.begin(), xs.begin(), xs.end(),
                  ys.begin(), ys.end(), xs_in_lcs);
    set_lcs(xs.begin(), xs_in_lcs, back_inserter(an_lcs));
}

    /// Update this element buffer into the GPU.
    ///
    /// Empty data will not be updated.
    ///
    /// \param data Data to update.
    void update(const seq<uint16_t> &data) const
    {
      if(!data)
      {
        return;
      }

#if 0
      for(unsigned ii = 0; (ii < data.size()); ii += 3)
      {
        std::cout << data[ii] << ", " << data[ii + 1] << ", " << data[ii + 2] << std::endl;
      }
#endif

      bind();
      dnload_glBufferData(GL_ELEMENT_ARRAY_BUFFER, data.getSizeBytes(), data.getData(), GL_STATIC_DRAW);
#if defined(USE_LD)
      vgl::increment_data_size_index(data.getSizeBytes());
#endif
    }

 bool operator>(const seq &b) const{
     if (ave > b.getAve()) {
         return true;
     }
     if (ave == b.getAve() && last>b.getLast()) {
         return true;
     }
     if (ave == b.getAve() && last == b.getLast() && end < b.getEnd()) {
         return true;
     }
     return false;
 }

   /// Compile this shader.
   ///
   /// Will terminate program on failure.
   void compile() const
   {
#if defined(USE_LD)
     std::string pretty_source = this->str();
     const GLchar *pretty_source_c_str = pretty_source.c_str();
     dnload_glShaderSource(m_id, 1, &pretty_source_c_str, NULL);
#else
     dnload_glShaderSource(m_id, m_parts.size(), const_cast<const GLchar**>(m_parts.getData()), NULL);
#endif

     dnload_glCompileShader(m_id);

#if defined(USE_LD)
     std::string log = GlslShaderSource::get_shader_info_log(m_id);
     if(0 < log.length())
     {
       std::cout << GlslShaderSource::string_add_line_numbers(pretty_source) << std::endl;
       std::cout << log << std::endl;
     }

     if(!GlslShaderSource::get_shader_compile_status(m_id))
     {
       teardown();
       exit(1);
     }
#endif
   }

    virtual
    void onExport(iwriter & writer)
    {
        writer.putUnsignedInt(_columns.size());
        writer.putBool(_useHighestStd);

        for (uint i=0; i!=_columns.size(); ++i)
            writer.putUnsignedInt(_columns[i]);
    }

    ZScore(Node * const parent, const bool useHighestStd=false) :
        Node(parent),
        _columns(parent->output().size()),
        _useHighestStd(useHighestStd)
    {
        _s1 = new double[_columns.size()];
        _s2 = new double[_columns.size()];

        for (uint i=0;i<_columns.size();++i)
            _columns[i] = i;
    }

    /// Get pass at index.
    ///
    /// Will create pass if not available.
    ///
    /// \param idx Pass index.
    /// \return Reference to pass.
    ObjectReferenceSeq& getPass(unsigned idx)
    {
      unsigned req_pass_count = idx + 1;

      if(m_objects.size() <= req_pass_count)
      {
        m_objects.resize(req_pass_count);
      }

      return m_objects[idx];
    }

    /// Get a fresh animation state.
    ///
    /// Data contained in the returned animation state is undefined, but not invalid.
    ///
    /// \return Animation state.
    AnimationState& newAnimationState()
    {
      unsigned ret = m_current_animation_state++;

      if(m_animation_states.size() <= ret)
      {
        m_animation_states.emplace_back();
      }

      return m_animation_states[ret];
    }

    /// Constructor.
    ///
    /// \param type GL_VERTEX_SHADER or GL_FRAGMENT_SHADER.
    /// \param part1 Shader part.
    Shader(GLenum type, const char *part1) :
      m_id(dnload_glCreateShader(type))
    {
      m_parts.push_back(part1);

      this->compile();
    }

    /// Interpolate animation frame from two other frames.
    ///
    /// \param src Source frame (earlier).
    /// \param dst Destination frame (later).
    /// \param current_time Animation time.
    void interpolateFrom(const AnimationFrame &lhs, const AnimationFrame &rhs, float current_time)
    {
      unsigned bone_count = lhs.getBoneCount();
#if defined(USE_LD)
      if(rhs.getBoneCount() != bone_count)
      {
        std::ostringstream sstr;
        sstr << "cannot interpolate between frames of size " << bone_count << " and " << rhs.getBoneCount();
        BOOST_THROW_EXCEPTION(std::runtime_error(sstr.str()));
      }
#endif

      //std::cout << "resizing bones or not? " << m_bones.size() << " vs " << bone_count << std::endl;

      m_time = current_time;
      m_bones.resize(bone_count);

      for(unsigned ii = 0; (bone_count > ii); ++ii)
      {
        const BoneState &ll = lhs.getBoneState(ii);
        const BoneState &rr = rhs.getBoneState(ii);
        float ltime = lhs.getTime();
        float rtime = rhs.getTime();
        float mix_time = (current_time - ltime) / (rtime - ltime);
        vec3 mix_pos = mix(ll.getPosition(), rr.getPosition(), mix_time);
        quat mix_rot = mix(ll.getRotation(), rr.getRotation(), mix_time);

        //std::cout << "mix time: " << mix_time << ": " << mix_pos << " ; " << mix_rot << std::endl;

        m_bones[ii] = BoneState(mix_pos, mix_rot);
      }
    }

    /// Initialize from data.
    ///
    /// \param data Frame data.
    /// \param bone_amount Amount of bone elements.
    void initFromData(const int16_t *data, unsigned frame_amount, float scale)
    {
      m_time = fixed_8_8_to_float(data[0]);

      //std::cout << m_time << std::endl;

#if defined(USE_LD)
      if((frame_amount % 7) != 0)
      {
        std::ostringstream sstr;
        sstr << "invalid frame amount: " << frame_amount;
        BOOST_THROW_EXCEPTION(std::runtime_error(sstr.str()));
      }
#endif

      for(unsigned ii = 1; ((frame_amount + 1) > ii); ii += 7)
      {
        vec3 pos(static_cast<float>(data[ii + 0]),
            static_cast<float>(data[ii + 1]),
            static_cast<float>(data[ii + 2]));
        quat rot(fixed_4_12_to_float(data[ii + 3]),
            fixed_4_12_to_float(data[ii + 4]),
            fixed_4_12_to_float(data[ii + 5]),
            fixed_4_12_to_float(data[ii + 6]));

        //std::cout << pos << " ; " << rot << std::endl;

        m_bones.emplace_back(pos * scale, rot);
      }
    }

 /// Draw shadow caps for this state.
 ///
 /// \param prg Program to use.
 /// \param pass Which pass to draw.
 void drawShadowCaps(const Program &prg, unsigned pass = 0) const
 {
   if(m_objects.size() > pass)
   {
     for(const ObjectReference &vv : m_objects[pass])
     {
       vv.drawShadowCaps(prg);
     }
   }
 }

    /// Duplicate from given frame.
    ///
    /// \param op Frame to duplicate.
    void duplicate(const AnimationFrame &op)
    {
      unsigned bone_count = op.getBoneCount();

      m_time = op.getTime();
      m_bones.resize(bone_count);

      for(unsigned ii = 0; (bone_count > ii); ++ii)
      {
        m_bones[ii] = op.getBoneState(ii);
      }
    }

    virtual void
    onFinish()
    {
        Feature & f = output();

        for (uint j=0; j!=_columns.size(); ++j)
            _s1[j] = _s2[j] = 0.0;

        for (ulong i=0; i<_cache.size(); i+=f.size())
        {
            for (uint j=0;j<_columns.size();++j)
            {
                const double v = _cache[i+_columns[j]];
                _s1[j] += v;
                _s2[j] += v*v;
            }
        }

        if (_useHighestStd) {
            double highest = -1.0;
            for (uint j=0; j<_columns.size(); ++j)
            {
                meanNstd(_cache.size()/f.size(),
                        _s1[j], _s2[j], _s1[j], _s2[j]);

                if (_s2[j] > highest)
                    highest = _s2[j];
            }

            for (ulong i=0; i < _cache.size(); i+= f.size())
            {
                for (uint j=0;j<f.size();++j)
                    f[j] = _cache[i+j];

                for (uint j=0;j<_columns.size();++j)
                    f[_columns[j]] = (f[_columns[j]] - _s1[j]) / highest;

                publish(f);
            }
        }
        else
        {
            for (uint j=0;j <_columns.size(); ++j)
                meanNstd(_cache.size()/f.size(),
                        _s1[j], _s2[j], _s1[j], _s2[j]);

            for (ulong i=0; i<_cache.size(); i+=f.size())
            {
                for (uint j=0;j<f.size();++j)
                    f[j] = _cache[i+j];

                for (uint j=0; j<_columns.size(); ++j)
                    f[_columns[j]] = (f[_columns[j]] - _s1[j]) / _s2[j];

                publish(f);
            }
        }

        _cache.clear();
    }

    /// Initializes the complete intro direction.
    ///
    /// \param data Data blob to initialize from.
    void initialize(const int16_t *data)
    {
      const int16_t* iter = data;

      for(;;)
      {
        m_scenes.push_back(new Scene(iter));

        iter = next_segment(next_segment(iter + 4));

        if(Spline::is_segment_end(iter))
        {
          break;
        }
      }
    }

    /// Resolve the scene.
    ///
    /// Writes camera and eye positions.
    ///
    /// \param stamp Timestamp (frames).
    /// \param out_id Scene id.
    /// \param out_cpos Camera position.
    /// \param out_epos Eye position.
    /// \param out_totaltime Total time in this scene.
    void resolveScene(unsigned stamp, SceneEnum &out_id, vec3 &out_cpos, vec3 &out_epos,
        unsigned &out_totaltime) const
    {
      for(unsigned ii = 0; (m_scenes.size() > ii); ++ii)
      {
        Scene* vv = m_scenes[ii];

        if(vv->getLength() > stamp)
        {
          out_id = vv->getId();
          out_cpos = vv->resolveCamera(stamp);
          out_epos = vv->resolveEye(stamp);

          // Calculate downwards to get total time spent in this scene.
          unsigned total_time = 0;
          for(;;)
          {
            if(0 >= ii)
            {
              break;
            }
            --ii;

            if(m_scenes[ii]->getId() == out_id)
            {
              total_time += m_scenes[ii]->getLength();
            }
          }
          out_totaltime = total_time + stamp;
          //std::cout << out_cpos << " ; " << out_epos << " ; " << out_totaltime << std::endl;
          return;
        }
        stamp -= vv->getLength();
      }

#if defined(USE_LD)
      BOOST_THROW_EXCEPTION(std::runtime_error("timestamp outside boundary"));
#else
      out_id = NONE;
      out_cpos = vec3(0.0f, 0.0f, 0.0f);
      out_epos = vec3(0.0f, 0.0f, 0.0f);
      out_totaltime = 0;
#endif
    }

 /// Get number of bones.
 ///
 /// \return Number of bones.
 unsigned getBoneCount() const
 {
   return m_bones.size();
 }

 /// Output to stream.
 ///
 /// \param ostr Output stream.
 /// \return Output stream.
 std::ostream& put(std::ostream &ostr) const
 {
   return ostr << "AnimationFrame(" << m_time << "): " << m_bones.size() << " bones";
 }

int main( int argc, char **argv )

{
  if (argc < 2) {
    cerr << "Usage: " << argv[0] << " scene.obj ..." << endl;
    exit(1);
  }

  glutInit(&argc, argv);
  glutInitDisplayMode( GLUT_DOUBLE | GLUT_RGBA | GLUT_DEPTH );

  int shift = (argc < 4 ? 0 : atoi(argv[3]) * 450);

  glutInitWindowSize( windowWidth, windowHeight );
  glutInitWindowPosition( 50 + shift, 50 );
  glutCreateWindow( "toon shading" );

  GLenum status = glewInit();
  if (status != GLEW_OK) {
    std::cerr << "Error: " << glewGetErrorString(status) << std::endl;
    return 1;
  }

  glutDisplayFunc( display );
  glutReshapeFunc( reshape );
  glutIdleFunc( idle );
  glutKeyboardFunc( keyPress );
  glutSpecialFunc( specialKeyPress );

  // Set up world objects

  for (int i=1; i<argc; i++)
    objs.add( new wfModel( argv[i] ) );

  if (objs.size() == 0) {
    cout << "no objects defined" << endl;
    exit(1);
  }

  // Find world centre and radius

  worldCentre = vec3(0,0,0);
  for (int i=0; i<objs.size(); i++)
    worldCentre = worldCentre + objs[i]->centre;
  worldCentre = (1/(float)objs.size()) * worldCentre;

  worldRadius = 0;
  for (int i=0; i<objs.size(); i++) {
    float radius = (objs[i]->centre - worldCentre).length() + objs[i]->radius;
    if (radius > worldRadius)
      worldRadius = radius;
  }

  // Point camera to the model

  const float initEyeDistance = 5.0;

  eyePosition = (initEyeDistance * worldRadius) * vec3(0.7,0.3,1).normalize() + worldCentre;
  fovy = 2 * atan2( 1, initEyeDistance );

  // Set up renderer

  axes = new Axes();

  renderer = new Renderer( windowWidth, windowHeight );

  // Go

  glutMainLoop();

  return 1;
}

 /// Create a new child bone.
 ///
 /// \param pos Position.
 /// \param rot Rotation.
 /// \return Reference to the child bone created.
 void addChild(Bone* op)
 {
   m_children.push_back(op);
   op->setParent(this);
 }