/*
 * @(#)ORonnModel.java  1.0 June 2010
 *
 * Copyright (c) 2010 Peter Troshin
 *
 *        BioJava development code
 *
 * This code may be freely distributed and modified under the
 * terms of the GNU Lesser General Public Licence.  This should
 * be distributed with the code.  If you do not have a copy,
 * see:
 *
 *      http://www.gnu.org/copyleft/lesser.html
 *
 * Copyright for this code is held jointly by the individual
 * authors.  These should be listed in @author doc comments.
 *
 * For more information on the BioJava project and its aims,
 * or to join the biojava-l mailing list, visit the home page
 * at:
 *
 *      http://www.biojava.org/
 *
 */
package org.biojava.nbio.ronn;

import org.biojava.nbio.ronn.ModelLoader.Model;
import org.biojava.nbio.ronn.ModelLoader.Threshold;
import org.slf4j.Logger;
import org.slf4j.LoggerFactory;

import java.io.File;
import java.io.FileNotFoundException;
import java.io.PrintWriter;



/**
 * Fully re-factored version of RONN model. Based on the code in C version of
 * RONN.
 *
 * @author Peter Troshin
 * @version 1.0
 * @since 3.0.2
 */
public final class ORonnModel {

        private static final Logger logger = LoggerFactory.getLogger(ORonnModel.class);

        /**
         * Order probability, corresponds to disorder as 1-order
         */
        private final float disorder_weight;

        private final static int AA_ALPHABET = 19;
        private final static int maxR = 110;
        //private final static float coef = 1.0f;
        /**
         * Holds encoded query sequence
         */
        private final short[] seqAA;
        /**
         * Holds query sequence
         */
        private final char[] query;

        private final Model model;

        /**
         * Disorder scores for all residues
         */
        private float[] scores = null;

        final float[] detect() {

        scores = new float[query.length];
        int sResidue;
        int dIndex;
        int r;
        float est, fOrder, pDisor, fDisor;
        final float[][] Z = new float[seqAA.length][ORonnModel.maxR];
        final int[] Q = new int[seqAA.length];
        final Threshold thold = new ModelLoader.Threshold(model.modelNum);

        /*
         * 19 looks like a size of the sliding window. So for any sequences
         * shorted than 19 AA the score will be NaN. Original RONN segfault in
         * such condition
         */
        for (sResidue = 0; sResidue <= query.length - ORonnModel.AA_ALPHABET; sResidue++) {
                est = 0.0f;

                for (dIndex = 0; dIndex < model.numOfDBAAseq; dIndex++) {
                final float[] rho = align(sResidue, dIndex);// search for the
                // maximum alignment between ith peptide from the
                // query and the dIndex-th database sequence
                est += model.W[dIndex] * Math.exp((rho[1] - rho[0]) / rho[0]);
                }

                fOrder = (float) (Math.exp(-0.5 * Math.pow(est - thold.mu0, 2.0)
                        / thold.sigma0) / (Math.sqrt(6.28) * thold.sigma0));

                fDisor = (float) (Math.exp(-0.5 * Math.pow(est - thold.mu1, 2.0)
                        / thold.sigma1) / (Math.sqrt(6.28) * thold.sigma1));

                pDisor = (float) (disorder_weight * fDisor / ((1.0 - disorder_weight)
                        * fOrder + disorder_weight * fDisor));
                for (r = sResidue; r < sResidue + ORonnModel.AA_ALPHABET; r++) {
                Z[r][Q[r]] = pDisor;
                Q[r]++;
                }
        }

        for (sResidue = 0; sResidue < query.length; sResidue++) {
                est = 0.0f;
                float[] zRow = Z[sResidue];
                int numOfIterations = Q[sResidue];
                for (r = 0; r < numOfIterations; r++) {
                est += zRow[r];
                }
                scores[sResidue] = est / numOfIterations;
        }
        return scores;
        }

        public void getScores(final File outfile) throws FileNotFoundException {
        final PrintWriter output = new PrintWriter(outfile);
        if (scores == null) {
                synchronized (this) {
                if (scores == null) {
                        detect();
                }
                }
        }
        for (int i = 0; i < scores.length; i++) {
                output.printf("%c\t%f\n", query[i], scores[i]);
        }
        output.close();
        }

        // sResidue query sequence index and dIndex database sequence index
        private final float[] align(final int sResidue, final int dIndex) {
        int dResidue, r;
        float maxScore = -1000000;
        float rho1 = 0;
        int maxIdx = 0;
        float rho0 = 0;
        short[] dbAARow = model.dbAA[dIndex];
        int numOfIterations = model.Length[dIndex] - ORonnModel.AA_ALPHABET;
        for (dResidue = 0; dResidue <= numOfIterations; dResidue++) {
                // go though the database sequence for maximised alignment
                rho1 = 0.0f;
                for (r = 0; r < ORonnModel.AA_ALPHABET; r++) {
                // go through the query sequence for one alignment
                rho1 += RonnConstraint.Blosum62[seqAA[sResidue + r]][dbAARow[dResidue
                        + r]];
                }
                if (rho1 > maxScore) {
                maxScore = rho1;
                maxIdx = dResidue;
                }
        }
        for (r = 0; r < ORonnModel.AA_ALPHABET; r++) {
                rho0 += RonnConstraint.Blosum62[dbAARow[maxIdx + r]][dbAARow[maxIdx
                        + r]];
        }
        return new float[] { rho0, maxScore };
        }

        public ORonnModel(final String sequence, final Model model,
                final float disorder) throws NumberFormatException {
        this.disorder_weight = disorder;
        this.model = model;
        query = sequence.toCharArray();
        seqAA = new short[query.length];
        assert model != null;
        assert model.numOfDBAAseq > 0;
        for (int sResidue = 0; sResidue < sequence.length(); sResidue++) {
                seqAA[sResidue] = RonnConstraint.INDEX[query[sResidue] - 'A'];
                if ((seqAA[sResidue] < 0) || (seqAA[sResidue] > 19)) {
                        logger.error("seqAA[sResidue]={}({})", seqAA[sResidue], query[sResidue]);
                        System.exit(1);
                }
        }
        }
}