lisyarus/psemek
1
0
Fork 1
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Neural animation: covariance matrix adaptation (wip)

Browse source
This commit is contained in:
Nikita Lisitsa 2021-01-16 17:18:32 +03:00
parent a11c5d7c39
commit 0a1d192327
1 changed files with 169 additions and 24 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

193
examples/animation_2d.cpp
View file

@ -22,6 +22,9 @@
#include <psemek/util/to_string.hpp>
#include <psemek/util/pretty_print.hpp>
#include <eigen3/Eigen/Core>
#include <eigen3/Eigen/Cholesky>
#include <fstream>
#include <iomanip>
@ -419,20 +422,18 @@ struct controller
// inputs should be equal to number of bones * 7
// outputs should be equal to number of joints
static constexpr std::size_t inputs = 12 * 7;
static constexpr std::size_t user_inputs = 0;
static constexpr std::size_t outputs = 11;
static constexpr std::size_t inputs = 24 * 7;
static constexpr std::size_t outputs = 23;
static constexpr std::size_t layer1 = outputs;
// static constexpr std::size_t layer2 = 4;
// static constexpr std::size_t layer3 = outputs;
// static constexpr std::size_t param_count = layer1 * (inputs + 1) + layer2 * (layer1 + 1) + layer3 * (layer2 + 1);
static constexpr std::size_t param_count = layer1 * (inputs + user_inputs + 1);
static constexpr std::size_t param_count = (inputs + 1) * outputs;
static constexpr float max_output = 20.f;
geom::matrix<float, outputs, inputs + user_inputs> m1;
geom::matrix<float, outputs, inputs> m1;
geom::vector<float, outputs> t1;
// geom::matrix<float, layer1, inputs> m1;
@ -442,8 +443,8 @@ struct controller
// geom::matrix<float, layer3, layer2> m3;
// geom::vector<float, layer3> t3;
// DIRTY HACKERY!!!!!
float * params() { return (float *)&m1; }
Eigen::VectorXf to_eigen() const;
void from_eigen(Eigen::VectorXf const & v);
int generation = 0;
@ -455,7 +456,7 @@ struct controller
void reset();
geom::vector<float, outputs> apply(system const & s, geom::vector<float, user_inputs> const & ui) const;
geom::vector<float, outputs> apply(system const & s) const;
void read(std::istream & is);
void write(std::ostream & os) const;
@ -506,6 +507,25 @@ void controller::randomize(RNG && rng, float amplitude)
// visit_v(t3);
}
Eigen::VectorXf controller::to_eigen() const
{
Eigen::VectorXf vec;
vec.setZero(param_count);
auto p = std::copy(m1.coords, m1.coords + inputs * outputs, vec.data());
std::copy(t1.coords, t1.coords + outputs, p);
return vec;
}
void controller::from_eigen(Eigen::VectorXf const & v)
{
if (v.size() != param_count)
throw std::runtime_error("Wrong data size");
std::copy(v.data(), v.data() + inputs * outputs, m1.coords);
std::copy(v.data() + inputs * outputs, v.data() + (inputs + 1) * outputs, t1.coords);
}
template <typename RNG>
void controller::mutate(RNG && rng, float amplitude)
{
@ -544,14 +564,14 @@ void controller::mutate(RNG && rng, float amplitude)
void controller::reset()
{}
geom::vector<float, controller::outputs> controller::apply(system const & sys, geom::vector<float, user_inputs> const & ui) const
geom::vector<float, controller::outputs> controller::apply(system const & sys) const
{
if (sys.bones.size() * 7 != inputs)
throw std::runtime_error(util::to_string("Wrong number of inputs: should be ", sys.bones.size() * 7));
if (sys.joints.size() != outputs)
throw std::runtime_error(util::to_string("Wrong number of outputs: should be ", sys.joints.size()));
geom::vector<float, inputs + user_inputs> input;
geom::vector<float, inputs> input;
for (std::size_t i = 0; i < inputs / 7; ++i)
{
input[7 * i + 0] = sys.bones[i].position[0] - sys.bones[0].position[0];
@ -562,8 +582,6 @@ geom::vector<float, controller::outputs> controller::apply(system const & sys, g
input[7 * i + 5] = sys.bones[i].velocity[1];
input[7 * i + 6] = sys.bones[i].angular_velocity;
}
for (std::size_t i = 0; i < user_inputs; ++i)
input[inputs + i] = ui[i];
auto r1 = activation(m1 * input + t1);
// auto r2 = activation(m2 * r1 + t2);
@ -673,7 +691,7 @@ struct animation_2d_app
} mode = mode::train;
std::vector<controller> population;
std::size_t const population_size = 256;
std::size_t const population_size = 1024;
std::size_t const max_train_frames = 10.f / physics.dt;
std::size_t const max_train_variations = 1;
float const position_variation_amplitude = 0.f;
@ -681,9 +699,12 @@ struct animation_2d_app
std::size_t train_iterations = 0;
std::size_t const max_train_iterations = 1024*8;
float const randomize_amplitude = 10.f;
float const initial_variance = 10.f;
static constexpr auto mutation_amplitude = [](float t){ return std::pow(10.f, 1.f + geom::lerp(0.f, -2.f, t)); };
Eigen::VectorXf mean;
Eigen::MatrixXf covariance;
float best_score = 0.f;
bool const warm_start = false;
bool const enable_testing = true;
@ -909,19 +930,18 @@ void animation_2d_app::reset_state(system & sys) const
sys.joints.push_back(joint{i, (i + 1) % (2 * n + 4), 1.f, -1.f, geom::interval<float>{geom::rad(0.f), geom::rad(60.f)}});
//*/
/*
// Worm
int n = 12;
int n = 5;
for (int i = 0; i < n; ++i)
sys.bones.push_back(bone{{i + 0.5f, 0.f}, {1.f, 0.f}, {0.f, 0.f}, 0.f, 1.f, 0.25f});
for (int i = 0; i + 1 < n; ++i)
sys.joints.push_back(joint{i, i + 1, 1.f, -1.f, geom::interval<float>{geom::rad(-30.f), geom::rad(30.f)}});
//*/
/*
// Caterpillar
int n = 5;
int n = 8;
for (int i = 0; i < n; ++i)
sys.bones.push_back(bone{{i + 0.5f, 0.5f}, {1.f, 0.f}, {0.f, 0.f}, 0.f, 1.f, 0.25f});
for (int i = 0; i < n; ++i)
@ -1077,7 +1097,7 @@ float animation_2d_app::eval_score(controller const & c, random::generator rng)
physics.advance(dt, sel, [&c, &energy, dt](system const & physics){
std::vector<float> torque(physics.joints.size(), 0.f);
auto ctrl = c.apply(physics, {});
auto ctrl = c.apply(physics);
for (std::size_t i = 0; i < ctrl.dimension; ++i)
torque[i] = ctrl[i];
energy += geom::length(ctrl) * dt;
@ -1147,6 +1167,8 @@ float animation_2d_app::eval_score(controller const & c, random::generator rng)
return score;
}
// Old evolutionary training implementation
void animation_2d_app::do_train()
{
if (population.empty())
@ -1165,7 +1187,7 @@ void animation_2d_app::do_train()
else
{
for (auto & c : population)
c.randomize(rng, randomize_amplitude);
c.randomize(rng, initial_variance);
}
}
else
@ -1194,11 +1216,14 @@ void animation_2d_app::do_train()
}
for (std::size_t i = preserved + cross_new; i < population.size(); ++i)
new_population[i].randomize(rng, randomize_amplitude);
new_population[i].randomize(rng, initial_variance);
population = std::move(new_population);
}
(void)&controller::to_eigen;
(void)&controller::from_eigen;
std::vector<std::pair<float, std::size_t>> scores(population.size());
std::vector<async::future<void>> futures;
for (std::size_t i = 0; i < population.size(); ++i)
@ -1222,6 +1247,126 @@ void animation_2d_app::do_train()
best_score = scores.front().first;
}
/*
// Covariance matrix adaptation
void animation_2d_app::do_train()
{
if (population.empty())
{
population.resize(population_size);
std::ifstream is(cache_location, std::ios::binary);
if (warm_start && is)
{
controller w;
w.read(is);
for (auto & c : population)
{
c = w;
c.mutate(rng, initial_variance);
}
}
else
{
for (auto & c : population)
c.randomize(rng, initial_variance);
}
mean.setZero(controller::param_count);
for (auto const & c : population)
{
Eigen::VectorXf v = c.to_eigen();
// mean += v;
}
mean /= static_cast<float>(population.size());
covariance.setZero(controller::param_count, controller::param_count);
for (auto const & c : population)
{
Eigen::VectorXf const v = c.to_eigen() - mean;
covariance += v * v.transpose();
}
covariance /= static_cast<float>(population.size() - 1);
}
else
{
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wmaybe-uninitialized"
auto ldlt = covariance.ldlt();
Eigen::MatrixXf const l = ldlt.matrixL();
Eigen::VectorXf const d = ldlt.vectorD();
random::normal_distribution<float> n;
Eigen::VectorXf v(controller::param_count);
for (std::size_t i = 0; i < population.size(); ++i)
{
for (std::size_t j = 0; j < controller::param_count; ++j)
v[j] = n(rng);
for (std::size_t j = 0; j < controller::param_count; ++j)
{
// if (d[j] < 0.f) throw std::runtime_error(util::to_string("Covariance matrix is not positive-semidefinite: ", d[j]));
v[j] *= std::sqrt(std::max(d[j], 0.f));
}
population[i].from_eigen(l * v + mean);
}
#pragma GCC diagnostic pop
}
std::vector<std::pair<float, std::size_t>> scores(population.size());
std::vector<async::future<void>> futures(population.size());
for (std::size_t i = 0; i < population.size(); ++i)
{
futures[i] = bg.dispatch([&, i, rng = rng]() mutable {
scores[i] = {eval_score(population[i], rng), i};
});
}
bg.wait_all(futures.begin(), futures.end()).get();
std::sort(scores.begin(), scores.end(), [](auto const & p1, auto const & p2){ return p1.first > p2.first; });
std::size_t const selected = population.size() / 4;
auto weight = [selected](std::size_t rank, float){
return static_cast<float>(selected - rank);
};
float sum_weights = 0.f;
for (std::size_t i = 0; i < selected; ++i)
sum_weights += weight(scores[i].second, scores[i].first);
mean.setZero(controller::param_count);
for (std::size_t i = 0; i < selected; ++i)
{
float w = weight(scores[i].second, scores[i].first) / sum_weights;
mean += w * population[scores[i].second].to_eigen();
}
Eigen::MatrixXf new_covariance;
new_covariance.setZero(controller::param_count, controller::param_count);
for (std::size_t i = 0; i < selected; ++i)
{
float w = weight(scores[i].second, scores[i].first) / sum_weights;
Eigen::VectorXf const v = population[scores[i].second].to_eigen() - mean;
new_covariance += w * v * v.transpose();
}
// new_covariance *= static_cast<float>(population.size() - 1) / static_cast<float>(population.size());
float const mu = 0.0001f;
covariance = covariance * (1.f - mu) + new_covariance * mu;
best_score = scores.front().first;
++train_iterations;
(void)lerp;
}
*/
void animation_2d_app::do_test()
{
std::optional<geom::point<float, 2>> m;
@ -1295,7 +1440,7 @@ void animation_2d_app::do_test()
{
physics.advance(frame_clock.restart().count(), sel, [this](system const & physics){
std::vector<float> torque(physics.joints.size());
auto ctrl = population[test_id].apply(physics, {});
auto ctrl = population[test_id].apply(physics);
for (std::size_t i = 0; i < ctrl.dimension; ++i)
torque[i] = ctrl[i];
return torque;