#include "libscl.h"
using namespace scl;
using namespace std;

int main(int argc, char** argp, char** envp)
{
  string progname(argp[0]);

  ifstream fin;
  ofstream fout;
  string filename;

  if (argc != 5) {
    cerr << "Usage: " << progname 
         << " snp_parmfile_name " << " parm_1 parm_2 parm_3" << '\n';
    error("Error, bad input line");
  }

  filename = argp[1];
  fin.open(filename.c_str());
  if (!fin) error("Error, " + progname + ", cannot open " + filename);

  cout << '\n';
  cout << progname << ": input file name is " << filename << '\n';

  REAL parm_1;
  if (!isREAL(argp[2],parm_1)) {
    error("Error, " + progname + ", input for parm_1 is not a number");
  }

  REAL parm_2;
  if (!isREAL(argp[3],parm_2)) {
    error("Error, " + progname + ", input for parm_2 is not a number");
  }

  REAL parm_3;
  if (!isREAL(argp[4],parm_3)) {
    error("Error, " + progname + ", input for parm_3 is not a number");
  }

  string line;
  string first_word;
  string remaining_words;
 
  vector<string> pflines;

  string test_word = "Index";

  while (getline(fin,line)) {
    remaining_words = line;
    first_word = cutstr(remaining_words);
    if(cutstr(first_word) == test_word) break;
  }

  intvec  snp_index;
  intvec  index;
  realmat theta;
  realmat sdev;
  vector<string> descriptor;

  INTEGER count = 0;
  while (getline(fin,line)) { 

    string remaining_words = line;
    string first_word  = cutstr(remaining_words);
    string second_word = cutstr(remaining_words);
    string third_word = cutstr(remaining_words);
    cutstr(remaining_words);
    string fourth_word = cutstr(remaining_words);

    test_word = "A(1,1)";

    INTEGER i;
    if (isINTEGER(first_word,i)) {
      if (fourth_word != test_word) {
        snp_index.push_back(i);
        ++count;
        index.push_back(count);
	REAL r;
	if (isREAL(second_word,r)) {
	  theta.push_back(r);
        }
	else {
	  error("Error, " + progname + ", second_word is not REAL");
        }
	if (isREAL(third_word,r)) {
	  sdev.push_back(r);
        }
	else {
	  error("Error, " + progname + ", third_word is not REAL");
        }
	descriptor.push_back(fourth_word + remaining_words);
      }
    }
  }
  fin.clear(); fin.close();

  index.push_back(++count);
  index.push_back(++count);
  index.push_back(++count);

  INTEGER l_parms = count;

  snp_index.push_back(0);
  snp_index.push_back(0);
  snp_index.push_back(0);

  theta.push_back(parm_1);
  theta.push_back(parm_2);
  theta.push_back(parm_3);

  sdev.push_back(fabs(parm_1)*0.1);
  sdev.push_back(fabs(parm_2)*0.1);
  sdev.push_back(fabs(parm_3)*0.1);

  descriptor.push_back("cmf_parms[1]");
  descriptor.push_back("cmf_parms[2]");
  descriptor.push_back("cmf_parms[3]");

  filename += ".pf";
  fout.open(filename.c_str());
  if (!fout) error("Error, " + progname + ", cannot open " + filename);

  cout << progname << ": output file name is " << filename << '\n';
  cout << '\n';

  fout << "PARAMETER START VALUES (required)" << '\n';

  for (INTEGER i=1; i<=l_parms; ++i) {
    fout << fmt('e',26,6,theta[i]) << "     1" << fmt('d',5,index[i])
         << fmt('d',5,snp_index[i]) << "  " << descriptor[i-1] << '\n';
  }

  fout << "PROPOSAL SCALING (required)" << '\n';

  for (INTEGER i=1; i<=l_parms; ++i) {
    fout << fmt('e',26,6,sdev[i]) << "      " << fmt('d',5,index[i])
         << fmt('d',5,snp_index[i]) << "  " << descriptor[i-1] << '\n';
  }

  fout << '\n';
  fout << '\n';
  fout << "Code for npbusr prior" << '\n';
  fout << '\n';

  fout << "    " << "l_parms = " << l_parms << ';' << '\n';
  fout << '\n';
  fout << "    " << "mean.resize(l_parms,1);" << '\n';
  fout << "    " << "sdev.resize(l_parms,1);" << '\n';
  fout << '\n';
  for (INTEGER i=1; i<=l_parms; ++i) {
    fout << "    " << "mean[" << fmt('d',3,i) << "] = " 
         << fmt('e',24,16,theta[i]) << ';' 
	 << "       // " << descriptor[i-1] <<'\n';
  }
  fout << '\n';
  for (INTEGER i=1; i<=l_parms; ++i) {
    fout << "    " << "sdev[" << fmt('d',3,i) << "] = " 
         << fmt('e',24,5,sdev[i]) << ';' 
	 << "       // " << descriptor[i-1] <<'\n';
  }

  fout.clear(); fout.close();
  
  return 0;
}
