2016-11-08

gensim の tfidf で正規化（normalize）に苦しんだ話

最近先輩に勧められて python の gensim というライブラリを使い始めたのですが、試しに tfidf やってみたらどうやって正規化してるのかわからなかったから調べたって話です。

かなり細かいことなのですが、同じことに苦しむ人がもしかしたらいるかもってことで記事にすることにしました。

gensim とは？

gensim は python で提供されている自然言語処理のライブラリで、tfidf や、LSI や LDA みたいなトピックモデル、はたまた word2vec なんかも手軽に計算できる便利なツールです。これ2008年からあるらしいんですけど知らなかった...これの存在知ってたら僕の卒論の実装もっと楽になった気がする(--;)

まあ過去の話はさておき、このライブラリを試してみるべくまずはtfidfの計算をしようとしたわけです。

2016-07-20

AGC001 参加記

気づけば前ブログを書いてから4ヶ月経ち、もう学生ではなくなっていた..

6月末くらいから再び競プロに参加できるようになり、先日記念すべきAGC001に参加しました！

agc001.contest.atcoder.jp

遂にAtCoderも世界展開ということで、第一回に参加できて嬉しく思います。結果だけ見ると二完で悲しい感じですが(--;)

以下は復習できた問題にだけ触れます。

ちなみにeditorialはこちら

A - BBQ Easy

A: BBQ Easy - AtCoder Grand Contest 001 | AtCoder

基本的に「一番大きな串に肉をできるだけ多く刺したいなら、対となる串もできるだけ大きくする」という発想で解きました。まあソートして順番に出力していけばOK。

実装

N = int(raw_input())
L = map(int,raw_input().split())
L.sort(reverse = True)
ans = 0
for i in xrange(N):
    ans += min(L[2*i],L[2*i+1])
print ans

B - Mysterious Light

B: Mysterious Light - AtCoder Grand Contest 001 | AtCoder

一応解けたものの、時間使いすぎて結局二完に終わってしまいました。シミュレーションするのは無理だろうってことで $N=2,3,..,11$ ぐらいまで三角形描いて $O(1)$ 解を考えたのですが（なぜかあると思い込んでた）、結局時間を使っただけという..。

ただ、例を何個も描いて行く過程で光の軌跡の長さの変化がユークリッドの互除法っぽくなってることに気づいてそれで解けました。想定解法が美しすぎて感動しましたね。

実装

N,X = map(int,raw_input().split())
X = min(X,N-X)
ans = 0
N -= X
while X > 0:
    ans += N/X * X *3
    N,X= X,N%X
print ans

C - Shorten Diameter

C: Shorten Diameter - AtCoder Grand Contest 001 | AtCoder

コンテスト中最後に考えていて結局解けなかった問題。解説聞くとああーって感じです。これは解きたかったなあ。コンテスト中は木の直径をとりあえず求めて、直径の長さを持つパスの端を削ってまた直径計算して..まで考えて、いやうまくいかんやんどうしようって言ってたら終了しました。

直径が $K$ なら木の中心から各ノードへの距離が $K/2$ になるってのは言われてみれば確かにそうですね。勉強になりました。

実装については $K$ が偶数のときは各ノードが中心であると仮定して、中心からの距離が $K/2$ 以下となるノードの数を数えてます。 $K$ が奇数のときはエッジが張られてる一組のノードをペアとして、このエッジを無視して両ノードからの距離が $(K-1)/2$ 以下となるノードの数を数えました。

実装

from sys import setrecursionlimit
setrecursionlimit(100000)
def rad_reach(b,h,dist,E,rad):
    cnt = 1
    if dist < rad:
        for nxt in E[h]:
            if nxt == b:
                continue
            cnt += rad_reach(h,nxt,dist+1,E,rad) 
    return cnt
def solve():
    N,K = map(int,raw_input().split())
    pairs = []
    E = [[] for i in xrange(N)]
    for i in xrange(N-1):
        a,b = map(lambda x:int(x)-1,raw_input().split())
        E[a].append(b)
        E[b].append(a)
        pairs.append((a,b))
    ans = N
    rad = K/2
    if K % 2 == 0:
        for c in xrange(N):
            ans = min(ans,N-rad_reach(-1,c,0,E,rad))
    else:
        for c1,c2 in pairs:
            ans = min(ans,N-rad_reach(c2,c1,0,E,rad)-rad_reach(c1,c2,0,E,rad))
    print ans
solve()

再帰で書いたほうがすっきりするんですけど、pythonで出すとギリギリTLEしました（pypyなら通る）。ところがこれを繰り返しで書くとpythonでもギリギリ通りました。前もこんなことあった気がします。

D - Arrays and Palindrome

D: Arrays and Palindrome - AtCoder Grand Contest 001 | AtCoder

コンテスト中読んだもののそもそも問題の意味がわからなくて、解説放送見てやっと理解しました。なんか数列の問題かと思ったら一筆書きの話になっててびっくりした。

まだ完全に解説が理解できてないので解けたら復習記事書こうと思います。

E - BBQ Hard

E: BBQ Hard - AtCoder Grand Contest 001 | AtCoder

コンテスト中読んで意味はわかったけどどうしていいかわからず放置しました。問題文読んだ中で解説聞いて一番感動した問題です。なのでこの問題だけ詳しく書こうと思います。

ある二つの串焼きセット $i,j$ を選んだとき、串焼きの刺し方は同じものを含む順列を考えて

$\frac{(A_i+B_i+A_j+B_j)!}{(A_i+A_j)!(B_i+B_j)!}$ （通り）

になる（前準備しとけば $O(1)$ で求められる）ってのはいいんですが、問題はこれを全部の組やると ${}_NC_2$ 組だけ計算することになるということ。これだと計算量が $O(N^2)$ になるので余裕で死にます。それで諦めたのですが..

解法としては、 ${}_{n+m}C_m = \frac{(n+m)!}{n!m!}$ が座標平面上で $(0,0)$ から $(n,m)$ まで左と下に進まずに至る経路の数であることを利用してます。解説放送見たとき「天才か！」って思いましたね。

サンプルを例に取ると、串焼きセット $3:(2,1)$ と串焼きセット $1:(1,1)$ で作れる串焼きの種類は $(0,0)$ から $(2+1,1+1)$ へ至る経路の数になります。で、解説だと分かりやすいように、 $(-2,-1)$ から $(1,1)$ への経路数と言い換えてます。

この経路をdpを使って求める方法は計算量が $O(1)$ になりませんが、重要なのは複数の組での計算が同時に行えることです。串焼きセット3と1,串焼きセット2と1で作れる串焼きの種類を以下のように一回のdpで一気に求められます。同じものを含む順列が出てくると条件反射で $O(1)$ で求めようとしちゃうので、これは本当にビックリしました！

f:id:Tawara:20160720003248p:plain f:id:Tawara:20160720003712p:plain f:id:Tawara:20160720003745p:plain

で、実際の実装としては全部の組み合わせを計算するために、各串焼きセット $(A_i,B_i)$ について $(-A_i,-B_i)$ に 1 を加算し、dpで計算した後は各 $(A_i,B_i)$ で値を回収して和を取ります。サンプルだと $(1,1)$ が二つあるので、 $(-1,-1)$ には2が加算され、 $(1,1)$ で値を二回回収します。

f:id:Tawara:20160720011914p:plain

すると求める値は $22 \times 2 + 35 = 79$ となる..あれ？本来の解 $10 + 10 + 6 = 26$ と全然違います。理由は下の表に示したいらない部分(青字と黒字)の値も含まれているからです。

f:id:Tawara:20160720011349p:plain

なので、同じ串焼きセットで串焼きを作った場合の種類数を引いて $2$ で割れば、 $(79 - (15+6+6)) \div 2 = 26$ と答えが出せます。これで今度こそOKです。

余談

ところでいきなりなのですが、以下に示すpythonコードの一部を見てください。( $p = 10^9+7$ です)

print ((ans - dup) / 2 + p) % p

違和感を覚えたでしょうか？僕はこれのせいで解法を理解したにもかかわらずめちゃくちゃ時間を無駄にしました。

当たり前なんですけど、modとってるのに割り算しちゃダメですよね..本当に何で気付かなかったんだろう..階乗の逆元は前準備のためにしっかり求めてるのに2はそのまま割るという愚行を犯しました。今度から気を付けます。

実装

pypyでも誤差レベルでギリギリ間に合ったのですが（1998ms）、今のPCにC++導入したのでC++でも解いてみました。※include等は省略

#define range(i,a,b) for(int i=(a); i < (b); i++)
#define rep(i,n) range(i,0,n)
 
#define BASE 2000
#define MAX 4000
#define MOD 1000000007
 
LL F[MAX*2+1],Finv[MAX+1];
LL dp[MAX+2][MAX+2];
int A[200000],B[200000];
 
LL pow_mod(LL b,LL e){
    LL ret = 1;
    for (;e>0;e>>=1){
        if(e&1) ret = ret * b % MOD;
        b = b * b % MOD;
    }
    return ret;
}
int main(){
    int N,A_MAX=0,B_MAX=0;
    LL ans = 0,dup;
    F[0] = 1;
    rep(i,MAX*2) F[i+1] = F[i]*(i+1) % MOD;
    Finv[MAX] = pow_mod(F[MAX],MOD-2);
    rep(i,MAX) Finv[MAX-i-1] = Finv[MAX-i]*(MAX-i) % MOD;
    rep(i,MAX+1) rep(j,MAX+1) dp[i][j] = 0;
    scanf("%d",&N);
    rep(i,N){
        scanf("%d %d",A+i,B+i);
        A_MAX = max(A_MAX,A[i]);
        B_MAX = max(B_MAX,B[i]);
        dp[BASE-B[i]][BASE-A[i]] += 1;
    }
    range(i,BASE-B_MAX,BASE+B_MAX+1){
        range(j,BASE-A_MAX,BASE+A_MAX+1){
            dp[i+1][j] = (dp[i+1][j]+dp[i][j])%MOD;
            dp[i][j+1] = (dp[i][j+1]+dp[i][j])%MOD;
        }
    }
    rep(i,N){
        dup = F[A[i]*2+B[i]*2]*Finv[A[i]*2]%MOD*Finv[B[i]*2]%MOD;
        ans = (ans + dp[BASE+B[i]][BASE+A[i]] - dup + MOD) % MOD;
    }
    printf("%lld\n",ans*pow_mod(2,MOD-2)%MOD);
    return 0;
}

F - Wide Swap

F: Wide Swap - AtCoder Grand Contest 001 | AtCoder

問題文すら読めていないので何も書けません。時間内に誰も解けなかったそうです。

まとめ

世界向け第一回ということもあってか、どの問題も学びがあるめっちゃいいコンテストでした。こんな問題とけるようになりたいし、思いつくようにもなりたいです。

次回がいつになるかわかんないですけど絶対に参加したい。

2016-03-08

SRM684 div2 参加記

ブログ書くのもう三ヶ月ぶりくらいですね。

とりあえずどうにか卒業出来ることが決まり、めっちゃ久々にSRM に参加しました。

せっかくなので久々に参加記も書くことにします。

Easy (250) - Istr

問題概要

アルファベットの文字列 $S$ に対して、文字列に含まれる各アルファベットの数の二乗の総和をその文字列のスコアとする。この文字列から文字を $k$ 個削除してスコアを最小にせよ。

制約

$1 \le |S| \le 50$ , $0 \le k \le |S|$ , 文字列に含まれるのは "a - z" のアルファベットのみ

解法

含まれる数の多いアルファベットから一個ずつ減らしていけばよい。制約が緩いので、一回削除するごとにソートしてました。

実装

class Istr:
    def count(self, s, k):
        cnt = collections.Counter(s)
        l = cnt.values()
        l.sort()
        for i in xrange(k):
            if l[-1] == 0: break
            l[-1] -= 1
            l.sort()
        return sum(i**2 for i in l)

$0 \le k \le |S|$ だから if 文の判定が実は要らないってあとから気付きました。まあ実害は無いですが。

Med (600) - DivFreed2

問題概要

以下の条件を満たす長さ $n$ の数列 $a_n$ は何通り作れるか。 $10^9+7$ で割った余りを答えよ。

$1 \le i \le n$ について $1 \le a_i \le k$
$1 \le i \le n-1$ について $a_i \le a_{i+1}$ または $a_i \bmod a_{i+1} \ != 0$

制約

$1 \le n \le 10, 1 \le k \le 10^5$

解法

$a_n$ を先頭から順番に決めていくとき、 $a_{i+1}$ は直前の $a_i$ にのみ依存します。もう明らかにDPの匂いがプンプンする。

末尾の値が t であるような問題文の条件を満たす長さ (i+1) の数列の数を dp[i+1][t]とすると、dp[i+1][t]は以下のいずれかを満たす s に対する dp[i][s] の総和になります。

$s \le t$
$s \bmod t \ != 0$

問題は各 t に対して s がこれらの条件を満たしているかを全てのsについて確かめると、計算量が $O(n \ k^2)$ になってTLEすること。この部分で工夫が必要です。

この問題を解決するために累積和を用います。一つ目の条件を満たすsに対するdp[i][s]の総和については、先にdp[i][s]の累積和を小さい方から取っておけば $O(1)$ で答えられます。

ポイントとなるのは二つ目の条件の方の求め方です。一つ目の条件を満たさないかつ二つ目の条件を満たすsに対するdp[i][s]の総和は、

( $s > t$ かつ $s \bmod t \ != 0$ を満たす sに対する dp[i][s]の総和) =

( $s > t$ を満たす sに対する dp[i][s]の総和) - ( $s > t$ かつ $s \bmod t \ = 0$ を満たす sに対する dp[i][s]の総和)

と考えてやることで簡単に求められます。

( $s > t$ を満たす sに対する dp[i][s]の総和)については、一つ目の条件のときと同様に先にdp[i][s]の累積和を大きい方から取っておけば $O(1)$ で答えられます。

( $s > t$ かつ $s \bmod t \ = 0$ を満たす sに対する dp[i][s]の総和)については、要するに t の倍数であるような s に対する dp[i][s]の総和です。tの倍数を順番に確かめるのは $O(\log_t k)$ で行えるので、t+1 から k まで全部確かめるよりも計算量が格段に減ります。

以上により計算量は大体 $O(n k \log k)$ となるのでTLEを脱せます。ちょっとpythonつらそうだったので C++で出しました。

実装

using namespace std;

#define P 1000000007

long long dp[10][100001],acum_u[100002],acum_l[100002];

class DivFreed2 {
    public:
    int count(int n, int k) {
        int i,j,cnt;
        
        long long comp_sum, ans = 0;
        for (i = 0; i < k+1; i++) dp[0][i] = 1;
        dp[0][0] = 0;
        for (i = 1; i < n; i++) for (j = 0; j < k+1; j++) dp[i][j] = 0;
        for (i = 1; i < n; i++){
            for (j = 0; j < k+1; j++) acum_u[j] = acum_l[j] = dp[i-1][j];
            acum_u[k+1] = acum_l[k+1] = 0;
            for (j = 0; j < k; j++){
                acum_l[j+1] = (acum_l[j+1] + acum_l[j]) % P;
                acum_u[k-j-1] = (acum_u[k-j-1] + acum_u[k-j]) % P;
            }
            for (j = 1; j < k+1; j++){
                comp_sum = 0;
                for (cnt=j*2;cnt < k+1;cnt+=j){
                    comp_sum = (comp_sum + dp[i-1][cnt]) % P;
                }
                dp[i][j] += (acum_l[j]%P + ((acum_u[j+1] - comp_sum) % P + P))%P;
                dp[i][j] %= P;
            }
        }
        for (j = 1; j < k+1; j++) ans = (ans + dp[n-1][j]) % P;
        return ans;
    }
};

accumがacumになってるところに英語力の無さを感じる..

Hard (900) - Autohamil

問題概要

0か1のみからなる文字列を読んで動作するある決定性有限オートマトンを考える。このオートマトンの状態は 0 から n-1 の計 n 個で、状態 $i$ のときに0を読むと状態 z0[i] に、1を読むとz1[i]に遷移する。

このとき、全ての状態にちょうど一回だけなるような文字列が存在するか？

制約

$1 \le n \le 50, |z_0| = |z_1| = n, 0 \le z_{0 i} , z_{1 i} \le n-1$

感想

解けなかったので感想だけ。問題概要があまり概要になってませんが、要はこのオートマトンの状態遷移図を有向グラフと見たときに、ハミルトン路が存在するか？という問題。今気付いたんですが、よくみたら問題名もAuto''hamil''でちゃんとハミルトンって主張してました笑

正直全然分からなかったのとMedバグらせて時間無かったのもあって手がつけられませんでした。また挑戦したいと思います。

まとめ

結果はdiv2 50位でレートも50ほど上がったので久々にしては頑張った方だと思いたい。ただしょうもないミスでバグらせるのはいい加減止めたいなあと思います。

レートが1163になり、やっと青くなるチャンスが巡ってきました。次こそはdiv1に上がれるようがんばります！

2015-11-25

CODE FESTIVAL2015 決勝 (当日編) (A~E)

code-festival-2015-final.contest.atcoder.jp

なんか今更感ありますが、当日の決勝参加記をば。

本番ではA~Eまでの5問しか解けなくて悔しかったです。復習は復習編で。

A - コード川柳 (2:17でAC)

与えられた文字列の長さが順に5・7・5になってるか確かめるだけ。なのですが焦りすぎて1WAを出す失態を犯しました。

S,T,U = raw_input().split()
print "valid" if len(S) == len(U) == 5 and len(T) == 7 else "invalid"

B - ダイスゲーム (16:44でAC)

解法聴いてそっかーってなったやつ。ちゃんと考えれば良かったんですが、 $1 \leq N \leq 256$ なので $O(N^2)$ でも行けるやろってことでDPしました。ただこれもバグる要素ないはずなのに変に単一リストでやろうとして（そして結局やらない）時間かかってしまった。勿体無い。

N = int(raw_input())
dice = [[0]*(6*N+1) for i in xrange(N)]
for i in xrange(6):
    dice[0][i+1] = 1
for i in xrange(N-1):
    for j in xrange(6*(i+1),-1,-1):
        for k in xrange(6):
            dice[i+1][j+k+1] += dice[i][j]
print max(xrange(N,6*N+1), key= lambda x: dice[N-1][x])

C - 寿司タワー (26:01でAC)

シミュレーションしました。タワーの順番を前から見ていって、シャリ・ネタまたはネタ・シャリの順番になっていたらスルー、シャリまたはネタが連続しているときは、既に分解しているパーツがあればそれを使い、なければ寿司を一個分解して使わなかった方のパーツをストックします。

N = int(raw_input())
S = map(int,list(raw_input()))
n = s = cnt = inc = 0
while inc < 2*N-1:
    if S[inc]^S[inc+1]:
        inc += 2
    else:
        if S[inc] == 1:
            if n > 0:
                n -= 1
            else: 
                cnt += 1
                s += 1
        else:
            if s > 0:
                s -= 1
            else:
                cnt += 1
                n += 1
        inc += 1
print cnt

if文多い。解説の解法は頭良いしすっきりしてて良いですよね。

D - 足ゲームII (102:31でAC)

問題見た直後にとりあえずいもす法で累積和取ったんですが、そのあとどうすれば良いか分からない..「絶対 $O(N)$ 解あるやろ！」と思ったんですが思いつかず、パーカー欲しさからセグメントツリーで殴ることを決断。

N-1点とるために必要な人数はあるボタンをスルーしたときに必要になる人数の最小値なので、セグツリーで時間の区間ごとの必要人数を管理し、各ボタン( $S_i$ から $T_i$ の間押し続けると得点がとれる)についてそのボタンをスルーしたときに必要になる人数、すなわち

区間 $0 \leq t \lt S$ で必要な人数の最大値
区間 $S \leq t \lt T$ で必要な人数の最大値 - 1
区間 $T \leq t \lt 10^5$ で必要な人数の最大値

の3つのうちの最大値を求めます。その最小値が答えです。計算量は $O(NlogN)$ なのでC++なら余裕で間に合います。

しかしながらバグらせたせいで75分も時間取られました。しかもバグじゃなくて単にinit()呼び出してないだけだったんでめっちゃ悔やまれる..

int n;
int tree[100000*4];
int s[100000];
int t[100000];
void init(int N){
    n = pow(2,int(ceil(log2(100000))))-1;
    rep(i,(n+1)*2) tree[i] = 0;
}
int update(int h){
    if (h < n) tree[h] = max(update(h*2+1),update(h*2+2));
    return tree[h];
}
int query(int a, int b, int h, int l, int r){
    if (r <= a || b <= l) return 0;
    if (a <= l && r <= b) return tree[h];
    return max(query(a,b,h*2+1,l,(l+r)/2), query(a,b,h*2+2,(l+r)/2,r));
}
int main(){
    int N,ans,tmp;
    cin >> N;
    init(100000);
    rep(i,N){
        cin >> s[i];cin >> t[i];
        s[i]--; t[i]--;
        tree[n+s[i]]++;
        tree[n+t[i]]--;
    }
    rep(i,100000) tree[n+i+1] += tree[n+i];
    update(0);
    ans = N;
    rep(i,N){
        tmp = 0;
        tmp = max(tmp,query(s[i],t[i],0,0,n+1)-1);
        tmp = max(tmp,query(0, s[i], 0, 0,n+1));
        tmp = max(tmp,query(t[i], 100000, 0, 0, n+1));
        ans = min(ans,tmp);
    }
    cout << ans << endl;
    return 0;
}

E - ショートコーディング(133:00でAC)

これに関しては順当に考えたというか、--は+になる、!!!は!になる、-のあとに!が来ると-が無効化されるってのを考えて行って、それをそのままシミュレーションした感じです。順番に処理しないとまずいかと思ってわざわざスタックを用意してごちゃごちゃやってます。無効化された場合はポップし、最後まで残ったものを添字順にならべてます。

S = raw_input()[::-1]
m = e = 0
mstk = []
estk = []
for i,s in enumerate(S):
    if s == "-":
        if m == 0: 
            m += 1
            mstk.append(i)
        else: 
            m -= 1
            mstk.pop()
    else:
        if m == 1:
            m -= 1
            mstk.pop()
        if e < 2:
            e += 1
            estk.append(i)
        else:
            e -= 1
            estk.pop()
order = [(j,False) for j in mstk]+[(j,True) for j in estk]
order.sort(key = lambda x: x[0])
print "".join("!" if v else "-" for i,v in order)[::-1]

なんかCでもこんなコード書いたなあと思いつつ提出。AC出てパーカー獲得確定したときすごく嬉しかったです。

解説見てこんな長々と書いて順に処理する必要はないと知ってちょっとショックでしたが、解けたので良しということで。

残り時間

F,G,Hを眺めて、Fが解けそうな気がしたので取り組んでたんですが解けず。結局そのまま時間切れで終了しました。

まとめ

セグメントツリー勉強してて本当に良かったと思いましたね(バグッたけど)、おかげでパーカー取れました。ただ、区間系来てよくわからなかったときにとりあえずせグツリー持ち出しちゃうのは良くない気がする..

6問目解きたかったんですが5問で終わったのは残念でした。ここで6問目を解ける実力を付けたい。