From 0a1d19232734499a7c4b64ca7e4e280164029c6b Mon Sep 17 00:00:00 2001 From: lisyarus Date: Sat, 16 Jan 2021 17:18:32 +0300 Subject: [PATCH] Neural animation: covariance matrix adaptation (wip) --- examples/animation_2d.cpp | 193 +++++++++++++++++++++++++++++++++----- 1 file changed, 169 insertions(+), 24 deletions(-) diff --git a/examples/animation_2d.cpp b/examples/animation_2d.cpp index a3571af7..d70e9576 100644 --- a/examples/animation_2d.cpp +++ b/examples/animation_2d.cpp @@ -22,6 +22,9 @@ #include #include +#include +#include + #include #include @@ -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 m1; + geom::matrix m1; geom::vector t1; // geom::matrix m1; @@ -442,8 +443,8 @@ struct controller // geom::matrix m3; // geom::vector 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 apply(system const & s, geom::vector const & ui) const; + geom::vector 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 void controller::mutate(RNG && rng, float amplitude) { @@ -544,14 +564,14 @@ void controller::mutate(RNG && rng, float amplitude) void controller::reset() {} -geom::vector controller::apply(system const & sys, geom::vector const & ui) const +geom::vector 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 input; + geom::vector 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 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 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{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{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 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> scores(population.size()); std::vector> 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(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(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 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> scores(population.size()); + std::vector> 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(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(population.size() - 1) / static_cast(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> m; @@ -1295,7 +1440,7 @@ void animation_2d_app::do_test() { physics.advance(frame_clock.restart().count(), sel, [this](system const & physics){ std::vector 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;