/* k-means */

options mprint;

%let indir  = ~/in;  * 入力ファイル格納ディレクトリ ;
%let outdir = ~/out; * 出力ファイル格納ディレクトリ ;

%let clus_num = 3;     * クラスタの数をここで指定する ;
%let maxiter  = 100;   * 最大反復回数をここで指定する ;
%let mchg_cnt = 9999;  * クラスタ割り当ての変化判定用マクロ変数 ;

* クラスタの重心格納用グローバルマクロ変数 ;
%macro global_var;
    %do i = 1 %to &clus_num.;
        %global gx&i.;
        %global gy&i.;
    %end;
%mend global_var;

%global_var;


%macro km;

    * ファイル読込 ;
    data km_data;
        infile "&indir./km_data.csv" dsd missover firstobs = 2;
        attrib x    length = 8 label = "x座標";
        attrib y    length = 8 label = "y座標";
        attrib clus length = 8 label = "クラスタ";
        input x y;
        clus = 0;
        n = _n_; * 読込時のデータの並びを保持するための変数 ;
    run;

    * 各オブザベーションに乱数を振る ;
    data km_g;
        set km_data;
        call streaminit(100);
        rn = rand("uniform");
    run;

    * 与えられた乱数でソート ;
    * 初期重心として利用される ;
    proc sort data = km_g;
        by rn;
    run;

    /* クラスタ割り当ての変化が無くなるまで
       または最大反復回数まで繰り返し処理   */
    %let iter = 1;
    %do %while(&mchg_cnt. ne 0 and &iter. <= &maxiter.);
        %mac_g;
        %chg_clus;
        %cal_g;

        %let iter = %eval(&iter. + 1);
    %end;

    * 読込時と同じ並びにソート ;
    proc sort data = km_data;
        by n;
    run;

    * 結果ファイル出力 ;
    ods csv file = "&outdir./km_result.csv";
        proc print data = km_data noobs;
            var x y clus;
        run;
    ods csv close;

%mend km;


%macro mac_g;

    * 得られた重心をマクロ変数に格納する ;
    data _null_;
        set km_g(obs = &clus_num.);
        %do i = 1 %to &clus_num.;
            if _n_ eq &i. then do;
                call symput("gx&i.", compress(put(x, best.)));
                call symput("gy&i.", compress(put(y, best.)));
            end;
        %end;
    run;

%mend mac_g;


%macro chg_clus;

    /* 新しいクラスタの重心までの距離を計算し
       一番近いクラスタに割り当てを変更する   */
    data km_data;
        set km_data;
        old_clus = clus;
        chg_cnt = 0;
        %do i = 1 %to &clus_num.;
            d&i. = sqrt((&&gx&i.. - x) ** 2 + (&&gy&i.. - y) ** 2);
        %end;

        d_chk = min(of d1 - d&clus_num.);
        clus = 0;
        %do j = 1 %to &clus_num.;
            if d_chk eq d&j. then clus = &j.;
        %end;

        * クラスタ割り当ての変化を判定 ;
        if old_clus ne clus then chg_cnt = 1;
    run;

    * クラスタ割り当ての変化総数を計算 ;
    proc summary data = km_data;
        var chg_cnt;
        output out = chg_chk(drop = _type_ _freq_) sum = ;
    run;

    * クラスタ割り当ての変化総数をマクロ変数に格納 ;
    data _null_;
        set chg_chk;
        call symput("mchg_cnt", compress(put(chg_cnt, best.)));
    run;

%mend chg_clus;


%macro cal_g;

    proc sort data = km_data;
        by clus;
    run;

    * 各クラスタの重心を計算 ;
    proc summary data = km_data;
        by clus;
        var x y;
        output out = km_g(drop = _type_ _freq_) mean = ;
    run;

%mend cal_g;


* メイン処理実行 ;
%km;

proc fastclus data = km_data maxclusters = 3 maxiter = 100 out = km_result; 
    var x y;
run;