nariさんの作品(ソースコード)

//------------------------------------------------------------------------------
/// @file
/// @author   ハル研究所プログラミングコンテスト実行委員会
///
/// @copyright  Copyright (c) 2017 HAL Laboratory, Inc.
/// @attention  このファイルの利用は、同梱のREADMEにある
///             利用条件に従ってください。
//------------------------------------------------------------------------------

// nari 解答ファイル
// スコア 14,710

#include "Answer.hpp"

#include <vector>
#include <array>
#include <queue>
#include <algorithm>
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <map>
#include <complex>
#include <numeric>

/// プロコン問題環境を表します。
namespace hpc {

// 便利マクロ
#define REP(i,n) for(int i=0;i<(int)(n);i++)
#define FOR(i,a,b) for(int i=(a);i<(int)(b);i++)
#define FORR(i,a,b) for(int i=(b)-1;i>=(int)(a);i--)

// 点を複素数で表現
using Point = std::complex<float>;
// hpc::Vector2 との相互変換
Vector2 PointToVector2(const Point &p){
    return Vector2(p.real(), p.imag());
}
Point Vector2ToPoint(const Vector2 &v){
    return Point(v.x, v.y);
}

// 各UFOへの家の割当を表現する型
using Assignment = std::deque<int>;
using Assignments = std::vector<Assignment>;
// UFOの動きを表現する型
using Movement = std::queue<int>;
using Movements = std::vector<Movement>;
// UFOの動きをシミュレートして各ターンの位置情報を保存した型
using Simulate = std::vector<Point>;
using Simulates = std::vector<Simulate>;
// シミュレートの結果発生するアクションを保存した型
// <<turn, ufo_id>, action>
using ActionDetail = std::pair<std::pair<int,int>,Action>;
using ActionDetails = std::vector<ActionDetail>;

class Solver2 {
public:
    // パラメータから情報を取得
    const float houseRad = Parameter::HouseRadius;
    const float officeRad = Parameter::OfficeRadius;
    const float smallUfoRad = Parameter::SmallUFORadius;
    const float largeUfoRad = Parameter::LargeUFORadius;
    const float smallUfoSpeed = Parameter::SmallUFOMaxSpeed;
    const float largeUfoSpeed = Parameter::LargeUFOMaxSpeed;

    // 家の数
    int houseNum;
    // 農場の位置
    Point officePos;
    // 家の位置
    std::vector<Point> housePos;

    // 現在のターン
    int turn;
    // 現在のステージ番号
    int stageCount;

    // 計算結果の各UFOのシミュレート情報
    Simulates answerMove;
    // シミュレート結果のアクション情報
    ActionDetails answerAction;
    // アクション情報のイテレータ
    int answerActionIterator;

    Solver2(){
        turn = 0;
        stageCount = -1;
    }
    void init(const Stage& stage){
        // 初期化
        turn = 0;
        stageCount++;
        // Stageクラスからデータを取得
        officePos = Vector2ToPoint(stage.office().pos());
        houseNum = stage.houses().count();
        housePos.resize(houseNum);
        REP(i,houseNum){
            housePos[i] = Vector2ToPoint(stage.houses()[i].pos());
        }
        // デバッグプリント
        fprintf(stderr, "%d : %d\n", stageCount, houseNum);
        // 計算
        auto answer = solve();
        answerMove = answer.first;
        answerAction = answer.second;
        answerActionIterator = 0;
    }
    void moveItems(Actions &actions){
        // シミュレート結果のアクションをそのまま追加する
        while(answerActionIterator < (int)answerAction.size()){
            if(answerAction[answerActionIterator].first.first <= turn){
                actions.add(answerAction[answerActionIterator].second);
                answerActionIterator++;
            }else{
                break;
            }
        }
    }
    void moveUFOs(TargetPositions &targets){
        // 各UFOのシミュレートどおりに動かす
        REP(u,10){
            if(turn < (int)answerMove[u].size()){
                targets.add(PointToVector2(answerMove[u][turn]));
            }else{
                targets.add(PointToVector2(answerMove[u].back()));
            }
        }
        // ターンを進める
        turn++;
    }

    // 計算の本体
    std::pair<Simulates, ActionDetails> solve(){
        // 目標のターン数を定める
        int targetTurn = 40 + (houseNum - 40) * 50 / 60;
        // 計算が成功するまで繰り返す
        while(true){
            // 大UFOの動きを目標ターンに収まる限界で決める
            Assignments lass(2);    // 大UFOの家割当
            Movements lmov(2);      // 大UFOの動き
            Simulates lsim(2);      // 大UFOの位置シミュレート
            std::vector<bool> available(houseNum, true);    // どの家が未配達かを保存する配列
            // 一度、大UFO(id=0)の動きを仮決めする
            auto tmp = largeDp(targetTurn, available, officePos);
            lass[0] = tmp.first;
            lmov[0] = tmp.second;
            lsim[0] = simulateMove(lmov[0], false);
            // 配達状況をリセット
            available.assign(houseNum, true);
            // 大UFO(id=0)の最終地点から遠い方向に動くように、大UFO(id=1)の動きを決める
            tmp = largeDp(targetTurn, available, lsim[0].back());
            lass[1] = tmp.first;
            lmov[1] = tmp.second;
            lsim[1] = simulateMove(lmov[1], false);
            // 大UFO(id=1)の最終地点から遠い方向に動くように、大UFO(id=0)の動きを再度決める
            tmp = largeDp(targetTurn, available, lsim[1].back());
            lass[0] = tmp.first;
            lmov[0] = tmp.second;
            lsim[0] = simulateMove(lmov[0], false);
            // 計算メモを初期化する
            tspMemo.clear();
            evalMemo.clear();
            // 小UFOの家割当を目標ターンに収まる限界で決める
            Assignments sass = smallAssignLessTurn(available, targetTurn, lsim);
            // 目標ターンに収まるように配達できない場合、目標ターンを5増やして再度計算する
            if(sass.size()==0){
                targetTurn += 5;
                fprintf(stderr, "-> %d\n", targetTurn);
                continue;
            }
            // 小UFOの動きを計算する
            Movements smov(8);
            REP(i,8){
                smov[i] = smallTsp(sass[i], lsim).first;
            }
            // 動きをシミュレートした物を答えとして返す
            return simulateAll(lmov, lsim, smov, false, true);
        }
    }

    // 最良スコアには使っていない関数
    // 小UFOの初期配置を決めてから、ビームサーチを用いて改善する方針だった
    std::pair<Simulates, ActionDetails> solve2(
        int targetTurn, int largeAssignParam, int smallAssignParam){
        // 大UFOの動きを決める
        Assignments lass(2);
        Movements lmov(2);
        Simulates lsim(2);
        std::vector<bool> available(houseNum, true);
        if(largeAssignParam == 0){
            auto tmp = largeDp(targetTurn, available, officePos);
            lass[0] = tmp.first;
            lmov[0] = tmp.second;
            lsim[0] = simulateMove(lmov[0], false);
            available.assign(houseNum, true);
            tmp = largeDp(targetTurn, available, lsim[0].back());
            lass[1] = tmp.first;
            lmov[1] = tmp.second;
            lsim[1] = simulateMove(lmov[1], false);
            tmp = largeDp(targetTurn, available, lsim[1].back());
            lass[0] = tmp.first;
            lmov[0] = tmp.second;
            lsim[0] = simulateMove(lmov[0], false);
        }
        // REP(i,2){
        //     Movement mov = lmov[i];
        //     fprintf(stderr, "%d : ", i);
        //     while(mov.size()){
        //         fprintf(stderr, "%d, ", mov.front());
        //         mov.pop();
        //     }
        //     fprintf(stderr, "\n");
        // }
        // 小UFOの初期分割を決める
        tspMemo.clear();
        evalMemo.clear();
        Assignments sass(8);
        if(smallAssignParam == 0){
            sass = smallAssignEqNum(available);
        }else if(smallAssignParam == 1){
            sass = smallAssignEqAngle(available, rand() / 100.f);
        }
        // 小UFOの分割を改善する
        if(smallAssignParam >= 0){
            // sass = improveSass(lsim, sass);
            sass = beamSearchSass(lsim, sass);
        }
        // 小UFOの動きを計算する
        Movements smov(8);
        REP(i,8){
            smov[i] = smallTsp(sass[i], lsim).first;
        }
        // REP(i,8){
        //     Movement mov = smov[i];
        //     fprintf(stderr, "%d : ", i);
        //     while(mov.size()){
        //         fprintf(stderr, "%d, ", mov.front());
        //         mov.pop();
        //     }
        //     fprintf(stderr, "\n");
        // }
        // answer
        return simulateAll(lmov, lsim, smov, false, true);
    }

    // UFOの動きをシミュレートする関数
    // 小UFOにも使えるようにしたが、小UFOのシミュレートには使わなかった
    Simulate simulateMove(Movement mov, bool isSmall){
        // 計算結果
        Simulate sim;
        // UFOの半径
        float rad = isSmall ? smallUfoRad : largeUfoRad;
        // UFOの速度
        float speed = isSmall ? smallUfoSpeed : largeUfoSpeed;
        // UFOの現在位置(初期値は農場の位置)
        Point pos = officePos;
        while(!mov.empty()){
            // 現在位置を保存する
            sim.push_back(pos);
            // Movement型の値に応じて移動先へ動く
            int target = mov.front();
            Point tpos = officePos;
            float trad = officeRad;
            if(target >= 0){
                tpos = housePos[target];
                trad = houseRad;
            }
            Point dir = tpos - pos;
            float len = std::abs(dir);
            if(len > speed){
                dir = dir / len * speed;
            }
            pos += dir;
            // 移動先に到達した場合、Movement型を更新する
            if(std::abs(pos-tpos) <= rad + trad){
                mov.pop();
            }
        }
        // 最終位置を保存する
        sim.push_back(pos);
        return sim;
    }

    // 動的計画法で大UFOの動きを決める関数
    // farPos に指定した位置から遠ざかるように決める
    std::pair<Assignment, Movement> largeDp(int targetTurn, std::vector<bool> &available, Point farPos){
        const int N = Parameter::MaxHouseCount;
        // 農場、家の各位置間の移動にかかるターン数を保存する配列
        // 0 ~ n-1 を家、n を農場とする
        static float cost[N+1][N+1];
        // 動的計画法に用いる配列
        // dp[position][delivered count] -> min turn
        static float dp[N+1][N+1];
        static int bef[N+1][N+1];

        int n = houseNum;
        // costを計算する
        REP(i,n+1){
            cost[i][i] = 0;
            REP(j,i){
                // 2点間の距離を大UFO速度で割り、切り上げたものをcostとする
                Point p1 = i==n ? officePos : housePos[i];
                Point p2 = j==n ? officePos : housePos[j];
                float dist = std::abs(p1-p2);
                cost[i][j] = cost[j][i] = ceil(dist / largeUfoSpeed);
            }
        }
        // 家を農場からの距離でソートする
        // 大UFOは農場から離れるようにのみ動くと仮定している
        std::vector<int> ids(n);
        REP(i,n)ids[i] = i;
        sort(ids.begin(), ids.end(), [&](int a, int b){
            return std::norm(housePos[a]-officePos) < std::norm(housePos[b]-officePos);
        });
        // dp配列の初期化
        REP(i,n+1)REP(j,n+1){
            dp[i][j] = 2521.f;
        }
        dp[0][0] = 0;
        // 動的計画法の計算
        REP(h,n){   // 次に訪れる家のIDでループ(農場からの距離でソートされている)
            int hid = ids[h];   // 家のID
            if(!available[hid]){
                continue;
            }
            REP(i,h+1)REP(j,n){
                // dp[i][j] から dp[h+1][j+1] を更新する
                int id = i==0 ? n : ids[i-1];
                float len = cost[hid][id];
                if(dp[i][j] + len < dp[h+1][j+1]){
                    dp[h+1][j+1] = dp[i][j] + len;
                    // 状態の前情報を記録しておく
                    bef[h+1][j+1] = i;
                }
            }
        }
        // 最大評価となる状態を取得する
        // 評価は 配達した家の数 - かかったターン数 / 目標ターン + 農場を中心としたfarPosと最終地点のなす角(rad)
        int acq = 0, ci = 0;
        float bestHyoka = -1e9;
        FOR(i,1,n+1)REP(j,n+1){
            if(dp[i][j] < targetTurn){
                float hyoka = j - dp[i][j] / targetTurn;
                if(std::norm(farPos-officePos) > 1.f){
                    Point v1 = farPos - officePos;
                    Point v2 = housePos[ids[i-1]] - officePos;
                    float angle = std::abs(std::arg(v2/v1));
                    hyoka += angle * 1.0f;
                }
                if(hyoka > bestHyoka){
                    acq = j; ci = i;
                    bestHyoka = hyoka;
                }
            }
        }
        // 最大評価となる状態への経路を取得する
        std::vector<int> path;
        while(ci > 0){
            int b = bef[ci][acq];
            int hid = ids[ci-1];
            path.push_back(hid);
            available[hid] = false;
            acq--;
            ci = b;
        }
        std::reverse(path.begin(), path.end());

        // 家割当と動きを計算する
        Assignment ass;
        Movement mov;
        for(int x : path){
            mov.push(x);
            ass.push_back(x);
            available[x] = false;
        }
        return std::make_pair(ass, mov);
    }

    // 最良スコアには使っていない関数
    // 家を農場中心で偏角ソートして等しい個数で小UFOに割り当てる
    Assignments smallAssignEqNum(std::vector<bool> available){
        int n = 0;
        std::vector<int> ids;
        REP(i,houseNum)if(available[i]){
            ids.push_back(i);
            n++;
        }
        std::sort(ids.begin(), ids.end(), [&](int a,int b){
            return std::arg((housePos[a]-officePos) / (housePos[b]-officePos)) < 0;
        });
        const int m = 8;
        Assignments ret(m);
        int iter = 0;
        REP(i,m){
            Assignment ass;
            int num = (n/m) + (i<n%m ? 1 : 0);
            while(num--){
                ass.push_back(ids[iter++]);
            }
            ret[i] = ass;
        }
        return ret;
    }

    // 最良スコアには使っていない関数
    // 家を農場中心で偏角ソートして等しい角度の範囲にある家を小UFOに割り当てる
    Assignments smallAssignEqAngle(std::vector<bool> available, float baseRad){
        int n = 0;
        std::vector<int> ids;
        REP(i,houseNum)if(available[i]){
            ids.push_back(i);
            n++;
        }
        std::sort(ids.begin(), ids.end(), [&](int a,int b){
            return std::arg(housePos[a]-officePos) < std::arg(housePos[b]-officePos);
        });
        const int m = 8;
        Assignments ret(m);
        const float pi = std::atan(1.f) * 4.f;
        int iter = 0;
        REP(i,m){
            float targetAngle = -pi + (i+1)*(2.f * pi)/m;
            if(i==m-1)targetAngle += 2521.f;
            Assignment ass;
            while(iter < (int)ids.size()){
                if(std::arg(housePos[ids[iter]]-officePos) < targetAngle){
                    ass.push_back(ids[iter++]);
                }else{
                    break;
                }
            }
            ret[i] = ass;
        }
        // 各割当が12件以下になるように調整する
        return rotateSass(ret, 12);
    }

    // 小UFOの割当を、目標ターン以下であるかぎり偏角順に増やしていく
    Assignments smallAssignLessTurn(std::vector<bool> &available, int targetTurn, Simulates &lsim){
        // 農場を中心に偏角ソート
        int n = 0;
        std::vector<int> ids;
        REP(i,houseNum)if(available[i]){
            ids.push_back(i);
            n++;
        }
        std::sort(ids.begin(), ids.end(), [&](int a,int b){
            return std::arg((housePos[a]-officePos) / (housePos[b]-officePos)) < 0;
        });
        std::reverse(ids.begin(), ids.end());   // pop_frontが無いためreverse
        const int m = 8;    // 小UFOの数
        Assignments ret(m);
        REP(i,m){
            Assignment ass;
            // 割当の数を二分探索する
            int low = 0, high = std::min(14, (int)ids.size())+1;
            while(low + 1 < high){
                int x = (low + high) / 2;
                Assignment tst;
                // 残っている家から偏角順に取っていく
                REP(j,x)tst.push_back(ids[ids.size()-1-j]);
                // 評価の結果目標ターンに収まるかどうかを判定
                if(evalTsp(tst, lsim) <= targetTurn){
                    low = x;
                }else{
                    high = x;
                }
            }
            // low が二分探索の結果
            // 割当を計算し、ids から削除する
            while(low--){
                ass.push_back(ids.back());
                ids.pop_back();
            }
            ret[i] = ass;
        }
        if(ids.size()){
            // 目標ターン以内で割り当てられなかったため空の割当を返す
            return Assignments();
        }
        return ret;
    }

    // 最良スコアには使っていない関数
    // 各小UFOの割当数が一定以下になるように割当を回転させる
    Assignments rotateSass(Assignments sass, int limit){
        int n = (int)sass.size();
        int cnt = 0;
        while(cnt < n){
            if((int)sass[cnt].size() <= limit){
                cnt++;
            }else{
                while((int)sass[cnt].size() > limit){
                    sass[(cnt+1)%n].push_front(sass[cnt].back());
                    sass[cnt].pop_back();
                }
                cnt = 0;
            }
        }
        return sass;
    }

    // ビームサーチに用いた状態
    struct State{
        float score, min, sum;  // 最終ターン数、小UFOの最小ターン数、各小UFOのターン数の総和
        int peek;               // ターン数が最大となっている小UFOのID、つまりボトルネック
        Assignments sass;       // 小UFOの割当
    };
    // 最良スコアには使っていない関数
    // 小UFOの割当と大UFOのシミュレート結果から State を計算する
    State getState(Assignments sass, Simulates lsim){
        State ret;
        ret.sass = sass;
        ret.score = 0.f;
        ret.sum = 0.f;
        ret.min = 1e9;
        ret.peek = 0;
        REP(i,(int)sass.size()){
            float x = evalTsp(sass[i], lsim);
            if(x > ret.score){
                ret.score = x;
                ret.peek = i;
            }
            ret.min = std::min(ret.min, x);
            ret.sum += x;
        }
        return ret;
    }
    // 最良スコアには使っていない関数
    // State の比較関数
    bool goodState(State &a, State &b){
        float asc = a.score * 8.f + a.sum;
        float bsc = b.score * 8.f + b.sum;
        if(asc != bsc) return asc < bsc;
        if(a.score != b.score)return a.score < b.score;
        if(a.sum != b.sum)return a.sum < b.sum;
        return a.min > b.min;
    }
    // 最良スコアには使っていない関数
    // 初期配置を受け取り、State に関してビームサーチを行う
    Assignments beamSearchSass(Simulates lsim, Assignments sass){
        const int BEAM_WIDTH = 2;
        const int BEAM_LENGTH = 20;
        std::vector<State> Q,R;
        Q.push_back(getState(sass, lsim));
        State bestState = Q[0];
        REP(_,BEAM_LENGTH){
            // ビーム幅の数だけ良い State を取得する
            int width = std::min(BEAM_WIDTH, (int)Q.size());
            std::partial_sort(Q.begin(), Q.begin()+width, Q.end(), [&](State a,State b){
                return goodState(a,b);
            });
            R.clear();
            // 各 State に対して近傍を計算する
            REP(i,width){
                State st = Q[i];
                if(goodState(st, bestState)){
                    bestState = st;
                }
                auto sas = st.sass;
                int id = st.peek;
                int n = (int)sas.size();
                // ボトルネックとなっている小UFOの割当を隣に移す
                for(int rest : sas[id])REP(to,n)if(to == (id+1)%n || to == (id+n-1)%n){
                    auto nxtSass = sas;
                    nxtSass[to].push_back(rest);
                    nxtSass[id].erase(std::find(nxtSass[id].begin(), nxtSass[id].end(), rest));
                    nxtSass = rotateSass(nxtSass, 14);
                    R.push_back(getState(nxtSass, lsim));
                }
            }
            std::swap(R,Q);
        }
        return bestState.sass;
    }

    // 最良スコアには使っていない関数
    // 小UFOの割当を少し変え、スコアがよくなるように貪欲に改善する(山登り法)
    Assignments improveSass(Simulates lsim, Assignments sass){
        int n = 8;
        Assignments bestSass = sass;
        std::vector<float> bestTotal;
        float bestScore;
        float bestMin;
        float bestSum;
        // init
        REP(i,n){
            bestTotal.push_back(approximateEval(bestSass[i], lsim));
        }
        bestScore = *std::max_element(bestTotal.begin(), bestTotal.end());
        bestMin = *std::min_element(bestTotal.begin(), bestTotal.end());
        bestSum = std::accumulate(bestTotal.begin(), bestTotal.end(), 0.f);
        REP(i,100){
            // modify
            int id = std::max_element(bestTotal.begin(), bestTotal.end()) - bestTotal.begin();
            if(sass[id].empty())continue;
            auto tmpsass = sass;
            REP(to,n)if(id!=to){
                for(int rest : tmpsass[id]){
                    sass = tmpsass;
                    sass[to].push_back(rest);
                    sass[id].erase(std::find(sass[id].begin(), sass[id].end(), rest));
                    sass = rotateSass(sass, 14);
                    // evaluate
                    std::vector<float> tot;
                    REP(j,n){
                        tot.push_back(approximateEval(sass[j], lsim));
                    }
                    float score = *max_element(tot.begin(), tot.end());
                    float min = *min_element(tot.begin(), tot.end());
                    float sum = std::accumulate(tot.begin(), tot.end(), 0.f);
                    if(score + sum < bestScore + bestSum || (score + sum == bestScore + bestSum && 
                        (sum < bestSum || (sum == bestSum && min > bestMin)))){
                        bestSass = sass;
                        bestTotal = tot;
                        bestScore = score;
                        bestMin = min;
                        bestSum = sum;
                    }
                }
            }
            sass = bestSass;
        }
        return bestSass;
    }

    // TSP型DPの最大の家の数
    const static int TSP_MAX = 14;
    // 計算結果のメモ
    std::map<Assignment, std::pair<Movement, float>> tspMemo;
    // 小UFOの割当と大UFOのシミュレート結果から、TSP型DPで最適な動きを計算する
    std::pair<Movement, float> smallTsp(Assignment ass, Simulates lsim){
        const int MAX_N = TSP_MAX;
        const int C = Parameter::SmallUFOCapacity;
        // 農場、家の各位置間の移動にかかるターン数を保存する配列
        // 0 ~ n-1 を家、n を農場とする
        static float cost[MAX_N+1][MAX_N+1];
        // TSP型DPのための配列
        // dp[k][j][i] := i番目の場所にいて、j個の荷物を持っていて、kの家に配達した場合の最小ターン数
        // 0 ~ n-1 を家、n を農場、n+1 を大UFO(id=0)、n+2 を大UFO(id=1) とする
        static float dp[1<<MAX_N][C+1][MAX_N+3];
        static int bef[1<<MAX_N][C+1][MAX_N+3];
        static int befj[1<<MAX_N][C+1][MAX_N+3];

        // メモされていたらその結果を返す
        if(tspMemo.count(ass)){
            return tspMemo[ass];
        }
        int n = ass.size();
        // calc を計算する
        REP(i,n+1){
            cost[i][i] = 0.f;
            REP(j,i){
                // 2点間の距離を小UFO速度で割り、切り上げたものをcostとする
                int id1 = i==n ? -1 : ass[i];
                int id2 = j==n ? -1 : ass[j];
                Point p1 = id1==-1 ? officePos : housePos[id1];
                Point p2 = id2==-1 ? officePos : housePos[id2];
                float dist = std::abs(p1-p2);
                cost[i][j] = cost[j][i] = ceil(dist / Parameter::SmallUFOMaxSpeed);
            }
        }
        // DP配列の初期化
        REP(k,1<<n)REP(j,C+1)REP(i,n+3){
            dp[k][j][i] = 2521;
        }
        dp[0][C][n] = 0;
        // DPの計算(もらうDP)
        FOR(k,1,1<<n){
            // 小UFOが家にいる場合(i = 0 ~ n-1)
            REP(i,n)if((k>>i)&1){
                int befMask = k ^ (1<<i);
                // 家から来る場合(from = 0 ~ n-1)
                REP(j,C){
                    befj[k][j][i] = j+1;
                    REP(from,n)if((k>>from)&1){
                        float nxtCost = dp[befMask][j+1][from] + cost[i][from];
                        if(nxtCost < dp[k][j][i]){
                            dp[k][j][i] = nxtCost;
                            bef[k][j][i] = from;
                        }
                    }
                }
                // 農場から来る場合(from = n, befj = C)
                {
                    int j = C-1, from = n;
                    befj[k][j][i] = j+1;
                    float nxtCost = dp[befMask][j+1][from] + cost[i][from];
                    if(nxtCost < dp[k][j][i]){
                        dp[k][j][i] = nxtCost;
                        bef[k][j][i] = from;
                    }
                }
                // 大UFOから来る場合(from = n+1, n+2, befj = C)
                FOR(from,n+1,n+3){
                    int j = C-1;
                    befj[k][j][i] = j+1;
                    float befCost = dp[befMask][j+1][from];
                    int befTurn = (int)ceil(befCost);
                    int target = from-n-1;
                    Point befPos = befTurn >= (int)lsim[target].size()
                                    ? lsim[target].back() : lsim[target][befTurn];
                    float nxtCost = befCost + ceil(std::abs(befPos - housePos[ass[i]]) / Parameter::SmallUFOMaxSpeed);
                    if(nxtCost < dp[k][j][i]){
                        dp[k][j][i] = nxtCost;
                        bef[k][j][i] = from;
                    }
                }
            }
            // 小UFOが農場にいる場合(i = n, j = C)
            {
                int i = n, j = C;
                // 家から来る場合(from = 0 ~ n-1)
                REP(fromj,C)REP(from,n)if((k>>from)&1){
                    float nxtCost = dp[k][fromj][from] + cost[i][from];
                    if(nxtCost < dp[k][j][i]){
                        dp[k][j][i] = nxtCost;
                        bef[k][j][i] = from;
                        befj[k][j][i] = fromj;
                    }
                }
            }
            // 小UFOが大UFOにいる場合(i = n+1, n+2, j = C)
            FOR(i,n+1,n+3){
                int target = i-n-1;
                int j = C;  // 荷物を最大限受け渡すことのみを想定する
                REP(fromj,2)REP(from,n)if((k>>from)&1){
                    float befCost = dp[k][fromj][from];
                    Point befPos = housePos[ass[from]];
                    // 大UFOの位置が不定のため、移動にかかるターン数を二分探索する
                    int low = 0, high = 32;
                    while(low+1<high){
                        int mid = (low+high)/2;
                        int trn = befCost + mid;
                        Point toPos = trn >= (int)lsim[target].size()
                                        ? lsim[target].back() : lsim[target][trn];
                        float len = std::norm(toPos - befPos) / smallUfoSpeed / smallUfoSpeed;
                        if(len <= mid * mid){
                            high = mid;
                        }else{
                            low = mid;
                        }
                    }
                    if(high == 32)continue;
                    float nxtCost = befCost + high;
                    if(nxtCost < dp[k][j][i]){
                        dp[k][j][i] = nxtCost;
                        bef[k][j][i] = from;
                        befj[k][j][i] = fromj;
                    }
                }
            }
        }
        // 最小ターンとなる状態を取得する
        int curi = 0, curj = 0;
        int curk = (1<<n)-1;
        REP(j,C+1)REP(i,n)if(dp[curk][j][i] < dp[curk][curj][curi]){
            curi = i;
            curj = j;
        }
        float total = dp[curk][curj][curi];
        // 最小ターンを達成する経路を取得する
        std::vector<int> path;
        while(curk>0){
            if(curi < n){
                path.push_back(ass[curi]);
            }else if(curi == n){
                path.push_back(-1);
            }else if(curi == n+1){
                // 大UFO(id=0) に動く場合は、1000 + ターン数 を保存する
                path.push_back(1000 + (int)dp[curk][curj][curi]);
            }else if(curi == n+2){
                // 大UFO(id=1) に動く場合は、2000 + ターン数 を保存する
                path.push_back(2000 + (int)dp[curk][curj][curi]);
            }
            int b = bef[curk][curj][curi];
            int bj = befj[curk][curj][curi];
            if(curi < n)curk ^= (1<<curi);
            curj = bj;
            curi = b;
        }
        std::reverse(path.begin(), path.end());
        // Movement に変換する
        Movement mov;
        for(int x : path){
            mov.push(x);
        }
        auto ret = std::make_pair(mov, total);
        // メモに保存しながら返す
        return tspMemo[ass] = ret;
    }

    // smallTsp 関数から Movement 計算を除き、ターン数計算のみにした関数
    std::map<Assignment, float> evalMemo;
    float evalTsp(Assignment ass, Simulates lsim){
        const int MAX_N = TSP_MAX;
        const int C = Parameter::SmallUFOCapacity;
        static float cost[32][32];
        static float dp[1<<MAX_N][16][32];

        if(evalMemo.count(ass)){
            return evalMemo[ass];
        }
        int n = ass.size();
        REP(i,n+1){
            cost[i][i] = 0.f;
            REP(j,i){
                int id1 = i==n ? -1 : ass[i];
                int id2 = j==n ? -1 : ass[j];
                Point p1 = id1==-1 ? officePos : housePos[id1];
                Point p2 = id2==-1 ? officePos : housePos[id2];
                float dist = std::abs(p1-p2);
                cost[i][j] = cost[j][i] = ceil(dist / Parameter::SmallUFOMaxSpeed);
            }
        }
        REP(k,1<<n)REP(j,C+1)REP(i,n+3){
            dp[k][j][i] = 2521;
        }
        dp[0][C][n] = 0;
        FOR(k,1,1<<n){
            REP(i,n)if((k>>i)&1){
                int befMask = k ^ (1<<i);
                REP(j,C){
                    float min = dp[k][j][i];
                    REP(from,n)if((k>>from)&1){
                        float nxtCost = dp[befMask][j+1][from] + cost[i][from];
                        min = std::min(min, nxtCost);
                    }
                    dp[k][j][i] = min;
                }
                {
                    int j = C-1, from = n;
                    float nxtCost = dp[befMask][j+1][from] + cost[i][from];
                    dp[k][j][i] = std::min(dp[k][j][i], nxtCost);
                }
                FOR(from,n+1,n+3){
                    int j = C-1;
                    float befCost = dp[befMask][j+1][from];
                    int befTurn = (int)ceil(befCost);
                    int target = from-n-1;
                    Point befPos = befTurn >= (int)lsim[target].size()
                                    ? lsim[target].back() : lsim[target][befTurn];
                    float nxtCost = befCost + ceil(std::abs(befPos - housePos[ass[i]]) / Parameter::SmallUFOMaxSpeed);
                    dp[k][j][i] = std::min(dp[k][j][i], nxtCost);
                }
            }
            {
                int i = n, j = C;
                float min = dp[k][j][i];
                REP(fromj,C)REP(from,n)if((k>>from)&1){
                    float nxtCost = dp[k][fromj][from] + cost[i][from];
                    min = std::min(min, nxtCost);
                }
                dp[k][j][i] = min;
            }
            FOR(i,n+1,n+3){
                int target = i-n-1;
                int j = C;
                REP(fromj,2)REP(from,n)if((k>>from)&1){
                    float befCost = dp[k][fromj][from];
                    Point befPos = housePos[ass[from]];
                    int low = 0, high = 32;
                    while(low+1<high){
                        int mid = (low+high)/2;
                        int trn = befCost + mid;
                        Point toPos = trn >= (int)lsim[target].size()
                                        ? lsim[target].back() : lsim[target][trn];
                        float len = std::norm(toPos - befPos) / smallUfoSpeed / smallUfoSpeed;
                        if(len <= mid * mid){
                            high = mid;
                        }else{
                            low = mid;
                        }
                    }
                    if(high==32)continue;
                    float nxtCost = befCost + high;
                    dp[k][j][i] = std::min(dp[k][j][i], nxtCost);
                }
            }
        }
        int curk = (1<<n)-1;
        float ans = 2521;
        REP(j,C+1)REP(i,n){
            ans = std::min(ans, dp[curk][j][i]);
        }
        return evalMemo[ass] = ans;
    }

    // 最良スコアには使っていない関数
    // DPを用いず、最近傍法によっておおよそのターン数を計算する
    float approximateEval(Assignment ass, Simulates lsim){
        int n = ass.size();
        Point pos = officePos;
        std::vector<bool> visited(n, false);
        float ans = 0.f;
        int items = n % 5;
        if(items==0)items = 5;
        REP(i,n){
            // 一番近い家を探す
            int to = -1;
            REP(j,n)if(!visited[j]){
                if(to==-1){
                    to = j;
                }else if(std::norm(housePos[ass[to]]-pos) > std::norm(housePos[ass[j]]-pos)){
                    to = j;
                }
            }
            // 移動する
            ans += ceil(std::abs(housePos[ass[to]]-pos) / smallUfoSpeed);
            pos = housePos[ass[to]];
            visited[to] = true;
            items--;
            if(i==n-1)break;
            if(items>0)continue;
            // 荷物が無くなったら農場か大UFOに動く
            float officeCost = ceil(std::abs(officePos-pos) / smallUfoSpeed);
            float l1cost = 1.f, l2cost = 1.f;
            Point l1pos, l2pos;
            while(true){
                int id = (int)l1cost < (int)lsim[0].size() ? (int)l1cost : (int)lsim[0].size()-1;
                l1pos = lsim[0][id];
                if(ceil(std::abs(l1pos-pos) / smallUfoSpeed) <= l1cost)break;
                l1cost += 1.f;
            }
            while(true){
                int id = (int)l2cost < (int)lsim[1].size() ? (int)l2cost : (int)lsim[1].size()-1;
                l2pos = lsim[1][id];
                if(ceil(std::abs(l2pos-pos) / smallUfoSpeed) <= l2cost)break;
                l2cost += 1.f;
            }
            if(officeCost <= l1cost && officeCost <= l2cost){
                ans += officeCost;
                pos = officePos;
            }else if(l1cost < l2cost){
                ans += l1cost;
                pos = l1pos;
            }else{
                ans += l2cost;
                pos = l2pos;
            }
            items = 5;
        }
        return ans;
    }

    // 最大速度制限がある状態で目的地に出来る限り近づくように移動し、その点を返す関数
    Point translate(Point p, Point to, float speed){
        Point dir = to - p;
        float len = std::abs(dir);
        if(len > speed){
            dir = dir / len * speed;
        }
        return p + dir;
    }

    // 大UFOの動きとシミュレート結果、小UFOの動きを受け取り、全体のシミュレート結果を返す関数
    // largeFreeがtrueの場合、lmov が空になった場合に近くの家に貪欲に配達しに行く
    // smallFreeがtrueの場合、smov が空になった場合に近くの家に貪欲に配達しに行く
    std::pair<Simulates, ActionDetails> simulateAll(Movements lmov, Simulates lsim, Movements smov, bool largeFree, bool smallFree){
        using pii = std::pair<int,int>;
        std::vector<bool> delivered(houseNum, false);   // 配達状況
        std::vector<Point> pos(10, officePos);          // 各UFOの現在位置 大UFO(id=0,1) 小UFO(id=2~9)
        std::vector<int> items(10, 0);                  // 各UFOの所持荷物
        int deliverCount = 0;       // 配達数
        int cturn = 0;              // 現在のターン数
        Simulates sims(10);         // シミュレート結果
        ActionDetails actions;      // アクション情報
        // 配達が終わるかターン数が2000未満の間繰り返す
        while(deliverCount < houseNum && cturn < 2000){
            // アクション
            REP(u, 10){
                bool isLarge = u<2;                                 // 大UFOか否か
                Movement &mov = isLarge ? lmov[u] : smov[u-2];      // UFOの動き
                float rad = isLarge ? largeUfoRad : smallUfoRad;    // 半径
                bool freee = isLarge ? largeFree : smallFree;       // 自由行動可能かどうか
                // 農場にいたらとりあえず積む
                bool officed = false;   // 農場で積んだかどうか
                if(std::abs(officePos - pos[u]) < rad + officeRad){
                    officed = true;
                    items[u] = isLarge ? 40 : 5;
                    actions.push_back(std::make_pair(pii(cturn, u), Action::PickUp(u)));
                }
                // Movement に従ったアクション
                if(!mov.empty()){
                    int pat = mov.front();
                    if(pat == -1){
                        // 農場に積みに行く
                        if(officed){
                            mov.pop();
                        }
                    }else if(pat >= 1000){
                        // 大UFOに積みに行く
                        int lu = (pat/1000) - 1;
                        // 積めるようなら積むアクションを作成
                        if(std::abs(pos[u]-pos[lu]) < rad + largeUfoRad){
                            actions.push_back(std::make_pair(pii(cturn, u), Action::Pass(lu, u)));
                            int passnum = std::min(5 - items[u], items[lu]);
                            items[u] += passnum;
                            items[lu] -= passnum;
                            mov.pop();
                        }
                    }else{  
                        if(freee && delivered[pat]){
                            // 自由行動可能で、配達済みならpop
                            mov.pop();
                        }else if(items[u] > 0 && std::abs(pos[u]-housePos[pat]) < rad + houseRad){
                            // そうでないなら配達する
                            if(!delivered[pat]){
                                actions.push_back(std::make_pair(pii(cturn, u), Action::Deliver(u, pat)));
                                items[u] --;
                                delivered[pat] = true;
                                deliverCount++;
                            }
                            mov.pop();
                        }
                    }
                }else{
                    // Movement が空なら重なっている家に配達する
                    REP(i,houseNum)if(!delivered[i]){
                        if(items[u] > 0 && std::abs(pos[u] - housePos[i]) < rad + houseRad){
                            actions.push_back(std::make_pair(pii(cturn, u), Action::Deliver(u, i)));
                            items[u] --;
                            delivered[i] = true;
                            deliverCount++;
                            break;
                        }
                    }
                }
            }
            // 移動
            REP(u, 10){
                bool isLarge = u<2;                                     // 大UFOか否か
                Movement &mov = isLarge ? lmov[u] : smov[u-2];          // UFOの動き
                bool freee = isLarge ? largeFree : smallFree;           // 自由行動可能かどうか
                float speed = isLarge ? largeUfoSpeed : smallUfoSpeed;  // UFOの速度
                // Movement に従った移動
                if(!mov.empty()){
                    int id = mov.front();
                    if(id == -1){
                        // 農場に向かう
                        pos[u] = translate(pos[u], officePos, speed);
                    }else if(id >= 1000){
                        // あるターンの大UFOの位置に向かう
                        int target = (id/1000) - 1;
                        int trn = id % 1000;
                        trn = std::min(trn, (int)lsim[target].size()-1);
                        if(std::abs(pos[u]-pos[target]) <= 40){     // 若干の位置ずれを考慮して実際の方向に動く
                            pos[u] = translate(pos[u], pos[target], speed);
                        }else{
                            pos[u] = translate(pos[u], lsim[target][trn], speed);
                        }
                    }else{
                        // 家に向かう
                        pos[u] = translate(pos[u], housePos[id], speed);
                    }
                }else if(freee){
                    // Movement が空で、自由行動可能ならば貪欲に配達をする
                    if(items[u] == 0){
                        // 荷物が無い場合は農場に向かう
                        pos[u] = translate(pos[u], officePos, speed);
                    }else{
                        // 荷物がある場合は一番近くの未配達の家に向かう
                        int to = -1;
                        REP(i,houseNum)if(!delivered[i]){
                            if(to == -1 || std::norm(pos[u]-housePos[to]) > std::norm(pos[u]-housePos[i])){
                                to = i;
                            }
                        }
                        if(to != -1){
                            pos[u] = translate(pos[u], housePos[to], speed);
                        }
                    }
                }else{
                    // 何もしない
                }
                sims[u].push_back(pos[u]);
            }
            cturn ++;
        }
        return std::make_pair(sims, actions);
    }
};

Solver2 g_Solver;
Answer::Answer(){}
Answer::~Answer(){}
void Answer::init(const Stage& aStage){g_Solver.init(aStage);}
void Answer::moveItems(const Stage& aStage, Actions& aActions){g_Solver.moveItems(aActions);}
void Answer::moveUFOs(const Stage& aStage, TargetPositions& aTargetPositions){g_Solver.moveUFOs(aTargetPositions);}
void Answer::finalize(const Stage& aStage){}

#undef REP
#undef FOR
#undef FORR

} // namespace
// EOF