自己相関（autocorrelation）

一般的に、場所や時間ごとの観測値の中に繰り返しパターンが潜んでいる場合、
「自己相関」を計算することでその周期性を検出できることがあります。
実際、音声やスペクトルの解析では、周期的に現れるシグナルを
ノイズの海の中から掬い上げるときに使われているそうです。

そして、自己相関は、化学構造にも適用されています。
ある距離（パス長, k）に位置する原子同士の性質（下のプログラムでは原子量）の積を
くまなく足しあわせ、値を算出しています。

　　Σ(wi x wj)
{i,j | d(i,j)=k }

一定の距離間隔で原子の性質に周期性が存在すれば、値が大きくなります。

k=0,1,2,3,4の5通りで、計算してみましょう。

/* autocor.java */
import java.io.*;

import org.openscience.cdk.interfaces.*;

import org.openscience.cdk.io.iterator.IteratingSMILESReader;

import org.openscience.cdk.interfaces.IMolecule;

import org.openscience.cdk.tools.manipulator.AtomContainerManipulator;
import org.openscience.cdk.config.IsotopeFactory;

import org.openscience.cdk.graph.matrix.TopologicalMatrix;

import org.openscience.cdk.interfaces.IAtomContainer;

import org.openscience.cdk.interfaces.IElement;

import org.openscience.cdk.tools.manipulator.AtomContainerManipulator;

import org.openscience.cdk.ChemObject;
import org.openscience.cdk.qsar.descriptors.molecular.AutocorrelationDescriptorMass;
class autocor {

    public static void main(String args[]){
        if(args.length!=1){

            System.err.println("autocor <smiles-file>");

            System.exit(1);

        }

        FileInputStream fis = null;

        IteratingSMILESReader isr = null;

        try{

            fis = new FileInputStream(new File(args[0]));

            isr = new IteratingSMILESReader(fis);

        } catch (Exception e) {

            e.printStackTrace();

        }
        while( isr.hasNext() ){

            IMolecule imol = (IMolecule)isr.next();

            String id = (String)imol.getProperty("cdk:Title");

            double[] masSum = new double[5];

            try {

                double[] w = listConvertion(imol);

                int natom = imol.getAtomCount();

                int[][] distancematrix = TopologicalMatrix.getMatrix(imol);

                for (int k = 0; k < 5; k++) {

                    for (int i = 0; i < natom; i++) {

                        for (int j = 0; j < natom; j++) {

                            if (distancematrix[i][j] == k) {

                                masSum[k] += w[i] * w[j];

                            } else masSum[k] += 0.0;

                        }

                    }

                    if (k > 0) masSum[k] = masSum[k] / 2;
                }

            } catch (Exception e1) {

                e1.printStackTrace();

            }
            System.out.printf("%s\t%.2f\t%.2f\t%.2f\t%.2f\t%.2f\n", id,

                              masSum[0], masSum[1], masSum[2], masSum[3], masSum[4]);

        }

    }


    private static double[] listConvertion(IAtomContainer container)

        throws java.io.IOException, ClassNotFoundException{

        double CARBON_MASS = 12.010735896788;

        int natom = container.getAtomCount();

        double[] scalated = new double[natom];

        for (int i = 0; i < natom; i++) {

            IsotopeFactory isofac = IsotopeFactory.getInstance(

                                           new ChemObject().getBuilder());

            String symbol = container.getAtom(i).getSymbol();

            scalated[i] = isofac.getNaturalMass( container.getAtom(i) ) / CARBON_MASS;

        }

        return scalated;

    }

}