ARC102 D - All Your Paths are Different Lengths

全然分かりませんでしたが、解説を見て理解できました。思いつけば簡単な問題なので悔しいです。

問題：D: All Your Paths are Different Lengths - AtCoder Regular Contest 102 | AtCoder

問題概要：パスの長さが0..L-1のちょうどL通りになるグラフを生成

解法： 例としてL=37の場合を考えてみます。

まず、長さが0から31までの32通りのパスからなるグラフは、以下のように構築できることに気付く必要があります。

L=37ですから、上記のグラフを改造して、長さが32, 33, 34, 35, 36のパスをさらに追加する必要があります。

ここで、頂点1から頂点3までは、ちょうど4種類の長さ0, 1, 2, 3が取れることに着目します。頂点3から頂点6までコスト32の辺をショートカットとして加えると、 ちょうど長さ32, 33, 34, 35のパスが追加できることがわかります。

残るは長さ36のパスです。今度は頂点1から頂点6までコスト36の辺を加えればよいです。よって以下が最終形です。

終了後ACしたコードを以下に記します。

#!/usr/bin/env python3
# -*- coding: utf-8 -*-

# 終了後に解説を見て解けた
# パス数が２のべき乗であれば簡単にグラフを構築できることに気付きたかった。
#
# 本番中は、10進法で考えて 頂点1と頂点2の間を辺{0, 100, 200, ..., 900}で繋ぎ、頂点2と頂点3の間を辺{0, 10, ..., 90}で繋ぎ、
# などという事を考えていたが、それ以上進展させることが難しかった。


import array
from bisect import *
from collections import *
import fractions
import heapq 
from itertools import *
import math
import random
import re
import string
import sys

L = int(input())

v = 1
edges = []  # (from, to, cost)

# L = 37の場合を例にとりコメントを付与
# まず長さ0 ~ 31までのパスを作成する。
#   頂点1 -> 頂点2を0の辺と1の辺で結ぶ
#   頂点2 -> 頂点3を0の辺と2の辺で結ぶ
#   頂点3 -> 頂点4を0の辺と4の辺で結ぶ
#   頂点4 -> 頂点5を0の辺と8の辺で結ぶ
#   頂点5 -> 頂点6を0の辺と16の辺で結ぶ
n = 2
cost = 1
while n <= L:
    edges.append((v, v+1, 0))
    edges.append((v, v+1, cost))
    cost *= 2
    v += 1
    n *= 2

# 頂点6をv_sink(最終到達点)とする
v_sink = v

# ここから先、長さ32, 33, 34, 35, 36のパスを生成することになる。
# 頂点3と頂点6を、32の辺で直でつなぐと、長さ32, 33, 34, 35のパスが作れる。
# 頂点1と頂点6を、36の辺で直でつなぐと、長さ36のパスが作れる。
n //= 2
N = n
L -= N

while L > 0:
    if L >= n:
        edges.append((v, v_sink, N))
        L -= n
        N += n
    n //= 2
    v -= 1

print(v_sink, len(edges))
for e in edges:
    print(*e)

SoundHound Programming Contest 2018 Masters Tournament 本戦 B - Neutralize

問題：B - Neutralize

問題概要：N個の薬品が並んでいる。各薬品の効用は-10⁹ ~ 10⁹。連続したK個の薬品の効用を0にする処理を任意回おこなったあとの薬品の効用の和の最大値を求める。

400点問題なのに2時間考えてまったくわかりませんでした。実行時間制限がなかったとしても解けてなかったと思います。sugim48さんの回答 ほぼそのまま参考にして本番後ACしたコードを最後に示します。また具体例をもとにアルゴリズムの動きを書き下した説明を以下に示します。

例としてK = 3で薬品が3 4 -8 2 -7 -5と並んだ場合を考えます。この解法では左側の薬品から順に部分配列を伸ばしていき、最後に増やした薬品を右端としたK個の薬品の効用を0にすべきかどうかを決めます。以下でその手順を見ていきましょう。

まず1個目を取り出して3。これではK個の薬品を取り出しようがないのでそのまま。この部分配列に対する答えは3。これをdp[1] = 3として記録します。

次に2個目を取り出して3 4。同様の理由でそのまま。この部分配列に対する答は7。これをdp[2] = 7として記録します。

次に3個目を取り出して3 4 -8。この部分配列に対する答は、「3 4 -8の3つの効用をゼロにするか、しないかのうち、良い方」です。今回はゼロにしたほうが得です。よってこの部分配列に対する答は0。これをdp[3] = 0として記録します。

次に4個目を取り出して3 4 -8 2。この部分配列に対する答は、「4 -8 2の3つの効用をゼロにするか、しないかのうち、良い方」です。4 -8 2の効用をゼロにしない場合は簡単で、dp[2]に2を足した値なので、2になります。4 -8 2の効用をゼロにする場合はちょっとややこしくて、4 -8 2だけでなくて4 -8 2のすぐ左側にある任意の数の薬品をもゼロにしてよい、と考えなければいけません。それを求めるにはこれまでのdpの値の最大値を記録する配列maを用意しておきます（詳細はコードを参照）。この場合は4 -8 2を0とするのがベストなので、3になります。2と3のうち良い方である3が答なので、dp[4] = 3として記録します。

次に5個目を取り出して3 4 -8 2 -7。この部分配列に対する答は、「-8 2 -7の3つの効用をゼロにするか、しないかのうち、良い方」です。上と同様の手順により、-8 2 -7をゼロにするのがベストなので、dp[5] = 7として記録します。

最後に6個目を取り出して3 4 -8 2 -7 -5。この部分配列に対する答は、「2 -7 -5の3つの効用をゼロにするか、しないかのうち、良い方」です。今回の場合は2 -7 -5の3つの効用をゼロにし、さらにその左側にある-8も巻き込んでゼロにするのがベストです。よってdp[6] = 7として記録します。

以上のアルゴリズムによるdp, maの値の変化を以下に示します。Bの値とdp, maの値との比較をしやすくするため、あえてdp, maの値は本来の位置よりひとつ左にずらしています。

idx	0	1	2	3	4	5
B	3	4	-8	2	-7	-5
dp	3	7	0	3	7	7
ma	3	7	7	7	7	7

しつこく考えてなんとか理解できたつもりですが、本当に正しいか、また人が読んでわかる説明になっているかは自信ありません。

CUDAを使って2枚の画像の平均をとるプログラム

CUDAの勉強がてら書いてみました。2枚の同サイズの画像を読み込んで、平均画像を作成して保存するだけのプログラムです。環境はUbuntu 16.04です。

// 2枚の画像の平均を取るプログラム
// 画像は以下から6000x4000のものを取得
//     https://www.pexels.com/photo/black-plants-photo-943903/
//     https://www.pexels.com/photo/person-riding-on-zip-line-1090551/

#include <chrono>
#include <cstdio>
#include <iostream>
#include <opencv2/core/core.hpp>
#include <opencv2/imgproc/imgproc.hpp>
#include <opencv2/highgui/highgui.hpp>

const int THREADS_PER_BLOCK = 16;

// https://developer.nvidia.com/cuda-example のサンプルから引用
static void HandleError( cudaError_t err, const char *file, int line ) {
    if (err != cudaSuccess) {
        printf( "%s in %s at line %d\n", cudaGetErrorString( err ), file, line );
        exit( EXIT_FAILURE );
    }
}
#define HANDLE_ERROR( err ) (HandleError( err, __FILE__, __LINE__ ))

__global__ void use_gpu_sub(unsigned char *a, unsigned char *b, unsigned char *c, int N) {
    int tid = threadIdx.x + blockIdx.x * blockDim.x;
    if (tid < N) {
        c[tid] = (a[tid] + b[tid]) / 2;
    }
}

void use_gpu(unsigned char* a, unsigned char* b, unsigned char* c, int height, int width, int channel) {
    unsigned char *dev_a, *dev_b, *dev_c;
    const int N = height * width * 3;

    // GPUのためのメモリ確保
    HANDLE_ERROR( cudaMalloc( (void**)&dev_a, N * sizeof(unsigned char) ) );
    HANDLE_ERROR( cudaMalloc( (void**)&dev_b, N * sizeof(unsigned char) ) );
    HANDLE_ERROR( cudaMalloc( (void**)&dev_c, N * sizeof(unsigned char) ) );

    // a, bがさす画像データをGPUのメモリにコピー
    HANDLE_ERROR( cudaMemcpy( dev_a, a, N * sizeof(unsigned char),
                            cudaMemcpyHostToDevice ) );
    HANDLE_ERROR( cudaMemcpy( dev_b, b, N * sizeof(unsigned char),
                            cudaMemcpyHostToDevice ) );

    // (N + THREADS_PER_BLOCK - 1)/THREADS_PER_BLOCK個のブロックをGPU上で実行。
    // 各ブロックでは THREADS_PER_BLOCK 個のスレッドが立つ
    use_gpu_sub<<<(N + THREADS_PER_BLOCK - 1)/THREADS_PER_BLOCK, THREADS_PER_BLOCK>>>( dev_a, dev_b, dev_c, N );

    // GPのメモリにある計算結果をCPU側にコピー
    HANDLE_ERROR( cudaMemcpy( c, dev_c, N * sizeof(unsigned char),
                            cudaMemcpyDeviceToHost ) );

    // free the memory allocated on the GPU
    HANDLE_ERROR( cudaFree( dev_a ) );
    HANDLE_ERROR( cudaFree( dev_b ) );
    HANDLE_ERROR( cudaFree( dev_c ) );
}

void use_cpu(unsigned char* a, unsigned char* b, unsigned char* c, int height, int width, int channel) {
    int N = height * width * channel;
    for (int i = 0; i < N; ++i) {
        c[i] = (a[i] + b[i]) / 2;
    }
}

void average_two_images() {
    cv::Mat img1 = cv::imread("adventure-asia-bali-1090551.jpg");
    cv::Mat img2 = cv::imread("beautiful-close-up-color-943903.jpg");

    assert(img1.size() == img2.size());
    assert(img1.type() == img2.type());
    std::cout << "image size: " << img1.size() << std::endl;

    cv::Mat img3_gpu = cv::Mat(img1.size(), img1.type());
    cv::Mat img3_cpu = cv::Mat(img1.size(), img1.type());

    // gpuを使って
    auto s1 = std::chrono::system_clock::now();
    use_gpu(img1.data, img2.data, img3_gpu.data, img1.rows, img1.cols, img1.channels());
    auto e1 = std::chrono::system_clock::now();
    std::cout << "gpu time: " << std::chrono::duration_cast<std::chrono::milliseconds>(e1 - s1).count() << " ms" << std::endl;

    // cpuを使って
    auto s2 = std::chrono::system_clock::now();
    use_cpu(img1.data, img2.data, img3_cpu.data, img1.rows, img1.cols, img1.channels());    
    auto e2 = std::chrono::system_clock::now();
    std::cout << "cpu time: " << std::chrono::duration_cast<std::chrono::milliseconds>(e2 - s2).count() << " ms" << std::endl;

    imwrite("averaged.jpg", img3_gpu);
    
    // 結果が同じことを確認
    assert(cv::countNonZero(img3_gpu != img3_cpu) == 0);
}

int main()
{
    average_two_images();
    return 0;
}

これを以下のコマンドでビルドします。

nvcc -O2 -std=c++11 average_two_images.cu $(pkg-config --cflags --libs opencv)

以下は実行結果です。GPUは転送速度がネックになっていてCPUより遅いです。

image size: [6000 x 4000]
gpu time: 180 ms
cpu time: 51 ms

コード中で定数にしている THREADS_PER_BLOCKの値を変えると処理時間が変わりました。どうやって最適な値を決めればよいのでしょう？

    176ms: 512スレッド
    172ms: 256スレッド
    158ms: 128スレッド
    158ms: 64スレッド
    180ms: 32スレッド
    180ms: 16スレッド

Visual Studio Codeでエディタ領域とターミナル領域とを行き来するショートカットを追加

Visual Studio CodeではCtrl + Shift + @というショートカットキーでターミナルを表示させることができます。ターミナル領域がアクティブな状態にてエディタ領域をアクティブにするにはCtrl + 1やCtrl + 2（数値はペインの番号を表す）が使えますが、逆に、エディタ領域がアクティブな状態にてターミナル領域をアクティブにするショートカットキーは、たぶんありません。強いて言えばCtrl + Shift + @を2度繰り返せばできますが、面倒です。

今のところのベストな方法が Switch focus between editor and integrated terminal in Visual Studio Code - Stack Overflow にありました。以下にその方法を転載します。

• Ctrl + Shift + P→Preferences: Open Keyboard shortcuts File を開く
• keybindings.jsonに以下の2行を追加

{ "key": "ctrl+`", "command": "workbench.action.terminal.focus"},
{ "key": "ctrl+`", "command": "workbench.action.focusActiveEditorGroup", "when": "terminalFocus"}

これにより、Ctrl + Shift + @を押すたびにフォーカスがエディタとターミナルを行ったり来たりするようにできました。

Docker imageのpushに失敗したらログインし忘れが原因かも

Dockerの練習中です。作成したdocker imageをdocker pushコマンドでDocker hubにpushしようとしたのですが、以下のエラーで失敗しました。

denied: requested access to the resource is denied

この場合、以下のことを確認する必要があります。

• https://hub.docker.com/ にアカウントを作成したか？
• docker loginコマンドでログインしたか？

自分の場合は後者が原因でした。

AnacondaのNumPyとPyPIのNumPyの速度を比較する

Anaconda Pythonで提供されるNumPyはIntelのMKLを利用しているため高速だという話を聞いたことがあります。実際どの程度違いがあるのか試してみました。

環境構築

実験は、自作PCに入れたUbuntu 16.04で行いました。環境構築にはAnaconda Pythonが提供する仮想環境を用います。

Anaconda版のNumPy

以下のコマンドで環境を作成しました。Cythonは必要ないかもしれませんが、chainerをインストールするときにmkl-random 1.0.1 requires cython, which is not installed. mkl-fft 1.0.2 requires cython, which is not installed.という警告が表示されたので、一応入れています。

$ conda create -n py37_conda_numpy python=3.7 numpy
$ source activate py37_conda_numpy
$ conda install cython
$ pip install chainer

インストールされたパッケージは以下のとおりです。

$ pip freeze
certifi==2018.4.16
chainer==4.3.0
Cython==0.28.3
filelock==3.0.4
mkl-fft==1.0.2
mkl-random==1.0.1
numpy==1.14.5
protobuf==3.6.0
six==1.11.0

PyPI版のNumPy

以下のコマンドで環境を作成しました。

$ conda create -n py37_pip_numpy python=3.7       
$ source activate py37_pip_numpy
$ pip install numpy
$ pip install chainer

インストールされたパッケージは以下のとおりです。

$ pip freeze
certifi==2018.4.16
chainer==4.3.0
filelock==3.0.4
numpy==1.14.5
protobuf==3.6.0
six==1.11.0

実験1. 行列積

1024x1024の行列二つの行列積を100回計算する以下のコードを使って、処理時間を比較しました。

import timeit
import numpy as np

def main():
    HEIGHT = 1024
    WIDTH = 1024
    REPEAT = 100

    np.random.seed(0)
    result = np.zeros((HEIGHT, WIDTH), dtype=np.float32)
    for _ in range(REPEAT):
        a = np.random.rand(HEIGHT, WIDTH)
        b = np.random.rand(HEIGHT, WIDTH)
        result += a @ b
    print(result.sum())

s = timeit.default_timer()
main()
print("elapsed time: {:.2f} sec.".format(timeit.default_timer() - s))

以下に結果を示します。

• Anaconda版NumPy: 3.35 秒
• PyPI版NumPy: 8.93 秒

Anaconda版のNumpyがPyPI版Pythonの38%程度の時間で処理を終えており、大差がつきました。プログラム実行中のCPUの負荷をhtopコマンドで眺めると、PyPI版のNumPyは8コアを全消費するのに比べてAnaconda版のNumPyは4コアしか消費しておらず、余裕を感じました。

実験2. ニューラルネット

手書きの数字を０から９に分類する分類器の学習chainer/train_mnist.py at master · chainer/chainer · GitHubを題材に、処理時間を比較しました。この分類器にはCNNは使われておらず、全結合層のみです。以下は実験の詳細です。

• ソースコードには実験1と同様の時間計測コードを追加した。
• 学習はCPUのみを用いる。GPUは用いない。
• 最初の実行時は実験用データセットをダウンロードする処理が走るので、それは処理時間に含めない。
• 引数に--epoch 5を指定し、5エポックで処理を打ち切る。

以下に結果を示します。

• Anaconda版NumPy: 54.51 秒
• PyPI版NumPy: 100.49 秒

Anaconda版のNumpyがPyPI版Pythonの54%程度の時間で処理を終えました。実験1ほどではないですが、やはりAnaconda版のNumpyは高速です。

まとめ

NumPyの速度がボトルネックになっているのであれば、Anaconda版のNumPyの導入は一考の価値があります。

AGC026 B rng_10s

おとといのAGCに参加して1問しか解けず。Bに二時間以上考えて解ききれないというのはさすがに問題だと思い、解説を読み込んで自分なりに咀嚼しました。本番で書いていたひどすぎるコードが消えてすっきりしました。しかし、本当に思考に穴がないかはまだ自信がありません。

こういう数学的な問題、いちど袋小路に入ってしまうと思考のループから抜け出せず時間を浪費してしまうことが多く、どうにもいけません。

問題：B - rng_10s 回答：ソースコード

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
#
# 本番中には解けなかった。解説pdfにある3つの場合分けはできた。
# B個の購入を可能なかぎりしたあとの個数が、「Cより大きくBより小さい」値に
# なることがあるとNoになる、ということには気づいたのだが、
# どう解いてよいか混乱して手がつけられなかった
#
#
# 解説pdf https://img.atcoder.jp/agc026/editorial.pdf を見た
#
# B個の購入を行った直後の個数の遷移を、mod Bの世界で考える。
#   最初は A mod B個。
#   B個の購入を可能な限り行った後も A mod B個で変わらず。
#   D個入荷すると    (A + D) mod B個。
# つまり、入荷するたびに、
# A からスタートして (A + D) mod B, (A + 2D) mod B, ... と遷移していく。
# 
# 遷移を無限に繰り返すと、mod Bの世界でのとりうる値は、g = gcd(B, D)として
#   (A % g), g + (A % g), 2g + (A % g), ..., kg + (A % g) (kは、kg + (A % g)がB未満となる最大の整数) のどれかとなる。
# ここでkを直接求めるのは難しい。
# かわりに、mを整数として、mg + (A % g)が最初にB以上となるのを求めるのは簡単で、それは B + (A % g)。
# これよりひとつ分小さいのは、B + (A % g) - gで、これが kg + (A % g)に等しい。
# 
# (A % g) + B - g が、C を超える場合、次に入荷が行われない。次もBを引くと個数が負になってしまうので No。
# (A % g) + B - g が、C 以下の場合、次に入荷が行われる。Yes。
#
#
# 解説PDFの、「また，(B/g − 1) × inv(D/g, B/g) 回 D を足したときにこの上界が達成できます」
# については理解できていない。
#



import array
from bisect import *
from collections import *
import fractions
import heapq 
from itertools import *
import math
import random
import re
import string
import sys


def gcd(a, b):
    while b:
        a, b = b, a % b
    return a


def solve(A, B, C, D):
    if B > A:
        return "No"
    if B > D:
        return "No"
    if B <= C:
        return "Yes"

    g = gcd(B, D)
    if B - g + (A % g) > C:
        return "No"
    else:
        return "Yes"


T = int(input())
for t in range(T):
    A, B, C, D = map(int, input().split())
    print(solve(A, B, C, D))