import java.io.*;
import java.util.*;
import org.biojava.bio.*;
import org.biojava.bio.seq.*;
import org.biojava.bio.seq.io.*;
import org.biojava.bio.seq.db.*;
import org.biojava.bio.dist.*;
import org.biojava.bio.dp.*;
import org.biojava.bio.symbol.*;
import org.biojava.bio.seq.Sequence;

public class PHMM 
{
  public PHMM()
  {
    try
    {
      /*
       * Make a profile HMM over the DNA Alphabet with 12 'columns' and default
       * DistributionFactories to construct the transition and emmission
       * Distributions
       */
      ProfileHMM hmm = new ProfileHMM(ProteinTools.getAlphabet(),
                                      12,
                                      DistributionFactory.DEFAULT,
                                      DistributionFactory.DEFAULT,
                                      "my profilehmm");

      //create the Dynamic Programming matrix for the model.
      DP dp = DPFactory.DEFAULT.createDP(hmm);
      
      //Database to hold the training set
      SequenceDB db = new HashSequenceDB();

      //Load the training set.....
      try
      { 
        //setup file path + name
        String dir             = "C:\\";
        BufferedInputStream is = new BufferedInputStream(new FileInputStream(dir + "input-fasta2.txt"));
        FileReader in          = new FileReader(dir + "input-proprietary.txt");        

        //get the appropriate alphabet:DNA,RNA, or PROTEIN
        Alphabet alpha = AlphabetManager.alphabetForName("PROTEIN");

        //get a sequenceDB of all sequences in the file
        db = SeqIOTools.readFasta(is,alpha);

      }catch (BioException ex) {
        //not in fasta format or wrong alphabet
        ex.printStackTrace();
      }catch (NoSuchElementException ex) {
        //no fasta sequences in the file
        ex.printStackTrace();
      }catch (FileNotFoundException ex) {
        //problem reading file
        ex.printStackTrace();
      }

      //train the model to have uniform parameters
      ModelTrainer mt = new SimpleModelTrainer();

      //register the model to train
      mt.registerModel(hmm);

      //as no other counts are being used the null weight will cause everything to be uniform
      mt.setNullModelWeight(1.0);
      mt.train();

      //create a BW trainer for the dp matrix generated from the HMM
      BaumWelchTrainer bwt = new BaumWelchTrainer(dp);

      //anonymous implementation of the stopping criteria interface to stop after 20 iterations
      StoppingCriteria stopper = new StoppingCriteria()
      {
        public boolean isTrainingComplete(TrainingAlgorithm ta)
        {
          return (ta.getCycle() > 20);
        }
      };
   
      /*
       * optimize the dp matrix to reflect the training set in db using a null model
       * weight of 1.0 and the Stopping criteria defined above.
       */
      bwt.train(db,1.0,stopper);

      /*
       * put the test sequence in an array, an array is used because for pairwise
       * alignments using an HMM there would need to be two SymbolLists in the 
       * array
       */
       
      //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      //THIS IS THE PART HAVING THE EXCEPTION THAT SAYS... 
      //"org.biojava.bio.symbol.IllegalSymbolException: Symbol d-1 not found in alphabet Transitions from d-1" 
      StatePath obs_rolls      = dp.generate(18);
      SymbolList roll_sequence = obs_rolls.symbolListForLabel(StatePath.SEQUENCE);
      SymbolList[] sla         = {roll_sequence};
      //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      
      //decode the most likely state path and produce an 'odds' score
      StatePath path = dp.viterbi(sla, ScoreType.PROBABILITY);
      System.out.println("Log Odds = "+path.getScore());

      //print state path
      for(int i = 1; i <= path.length(); i++)
      {
        System.out.println("g="+path.symbolAt(StatePath.STATES, i).getName());
      }
    }catch (Exception e){e.printStackTrace();}
  }

  public static void main(String[] args)
  {
    PHMM phmm = new PHMM();
    System.out.println("exit");
  }
}