// alicomp.cpp - Two alignments compare main prog - 29.05.2001 
//              


// Image size variables :
#define    _IMAGE_DEFAULT_WIDTH			640
//#define    _IMAGE_DEFAULT_WIDTH			800
#define    _IMAGE_DEFAULT_HEIGHT		480
#define    _IMAGE_MIN_WIDTH			    320
//#define    _IMAGE_MIN_HEIGHT		    240
//+#define    _IMAGE_MIN_HEIGHT		    150
#define    _IMAGE_MIN_HEIGHT		    200
#define    _IMAGE_MAX_WIDTH			    3000
//#define    _IMAGE_MAX_HEIGHT		    1500
#define    _IMAGE_MAX_HEIGHT		    10000

//#define    _IMAGE_DEFAULT_HEIGHT_DY		16
//#define    _IMAGE_MIN_HEIGHT_DY		    4
#define    _IMAGE_DEFAULT_HEIGHT_DY		20
#define    _IMAGE_MIN_HEIGHT_DY		    5

// ------------------------- a few debug flags
#define  _GETCH  //getch()
//#define  _TEST 
//#define  _TRACE 
//#define  _TRACE_ALIGN
//#define  _INTERNAL_NODES



// ------------------------------------- MSVC vs G++
#ifdef    _MSC_VER
#define   _MICROSOFT_
//+#include <config.h>    // included in calling .cpp

// Linux's g++
#else
#define CTRFLAG_USE_SGI_STL
#define __UNIX__
#define __ON_PC__
#define __LINUX__
// ? #define  _LINUX_
#endif

// --------------------- Some constants

#define   _ALIGN_BUFFER_LENGTH      1028
#define     _OPTIONS_STR_LEN            101
/**
#define   EZERO			0
#define   SLASH			PATH_SEPARATOR

#ifndef __ON_PC__
	#define M_SQRT2     1.41421356237309504880
#endif

#define	  STR_LEN		81
#define	  GAP			'-'
#define   GAP_STR		"-"

#define   A_HOMOLOGY 	16
#define   A_BACKGROUND 	32
#define   A_COLOR 	 	15
**/


// for graphic:
#define   _PI_			3.14159

// ----------------------------------------------

// Alex

// for cout
#include <iostream>

// for getch()
#ifdef _MSC_VER
	#include  <conio.h>
#endif

#include <stdio.h>
#include <stdlib.h>
#include <ctype.h>
#include <string.h>

// STL
#include <string>
//#include <ntree.h>
#include <vector>
#include <map>
//#include <multimap.h>
//#include <multiset.h>
#include <set>
using namespace std;

// WM from dhm.cpp  to compute column weights
//
// ------------------- __NEW logic (like swap logic) added to <javac.h>
//                       swap - to get around that guy in <javac.h>
#define   __NEW
#define   swap
/**
#include  "inv_sqrt.h"

#include  "const.h"
#include  "work_seq.h"
#include  "dothelix.h"
#include  "seqbase.h"
#include  "w_matrix.h"

#include  "seqbase.cpp"
#include  "inv_sqrt.cpp"
#include  "motif.cpp"
**/
#include  "seqbase.cpp"
#include  "params.cpp"
// ---------------------

//#include "w_matrix_pos.h"



typedef struct off_de {
	char	    ty;       // type: g - gap, 0(zero) - no delta, d - delta, e - end
                          //  ??   R - moved right, L - moved left (instead of delta) - ***NO NEED 

	int         ddx;      // delta; (is a Key abs(ddx) for the 2nd map )
    int         x1;       // starting from 0  (is a Key for the 1st map)
    int         x2;       // end pos

 // int         dir;      // 1 - +, -1 -

	off_de() : ty(' '), ddx(0) {}
	off_de( char c, int d = 0 ) : ty(c), ddx(d)
			{}
}   OffsetDescr;

// Alex
/*
typedef   vector< int, alloc >                                          vIntDef;
typedef   vector< string, alloc >                                       vAliDef;
typedef   vector< string, alloc >                                       vSNamesDef;
typedef   vector< map< int, OffsetDescr, less< int >, alloc >, alloc >  vOffsDef;
typedef   map< int, off_de, less< int >, alloc >                        mOffDef;
typedef   map< string, int, less< string >, alloc >                     mSNamesDef;
typedef   vector< int, alloc >                                          vRefDef;
*/
typedef   vector< int >                                          vIntDef;
typedef   vector< string >                                       vAliDef;
typedef   vector< string >                                       vSNamesDef;
typedef   vector< map< int, OffsetDescr, less< int > > >  vOffsDef;
typedef   map< int, off_de, less< int > >                        mOffDef;
typedef   map< string, int, less< string > >                     mSNamesDef;
typedef   vector< int >                                          vRefDef;

// used in aligfun.h
//typedef   vector< string, alloc >                                       vGroupDef;
typedef   vector< string >                                       vGroupDef;

// to call AligraImage() for type = 22
/*
typedef   vector< int, alloc >                                          vColGroupDef;
typedef   vector< double, alloc >                                       vGroupValueDef;
typedef   vector< int, alloc >                                          vSelGroupDef;
*/
typedef   vector< int >                                          vColGroupDef;
typedef   vector< double >                                       vGroupValueDef;
typedef   vector< int >                                          vSelGroupDef;


//static   char  *image_tail[] = { "a", "m", "se", "3", "4", "5", "6", "7", "8", "9", "10",
//                                 "s1", "s2"
//                                };  
typedef struct prog_options {
	//string  image_pref;            
	int     iW       ;             // width. Default 640 
	int     iH       ;             // height. Defaulu 480
// Alex 
//    vector< int, alloc >     vType;     // 0 - std, 1 - max in col, 2 - selected gr, 3 - most gr, 4 - least gr
    vector< int >     vType;     // 0 - std, 1 - max in col, 2 - selected gr, 3 - most gr, 4 - least gr
    //                                    // 11 - stat,   21 - comp (%%), 22 - comp 
    //vSelGroupDef             vSel;
    int     offType  ;    // offset scale type; 0 - from 0/1 (abs), 1 - from min (relative)

    string  image_pref;
    string  title;
    string  algo[2];
    string  sMatrixName;
    string  sMismatchFile;
    int     offMisMode;     // 0 - not display offsets for mismatched seq, 1 - display

    string  sGroupFile;    // color group file name - only for type = 22 (comp col values)
    string  sErrFile;      // molType err, ...
    int     fillType ;     // 0 - full (unweighted), 6 - weighted middle, 7 - very wei middle

    //                                             0 - AUTO adjust  
    prog_options( ) : iW(_IMAGE_DEFAULT_WIDTH), iH(0), title(""), //Graphical Image of Alignment"),
                      image_pref( "gali" ), //sMatrixName("blosum62.matrix"),
                      sMismatchFile("mismatch.txt"), offType(0), offMisMode(0),
                      //sGroupFile("color.def")
                      sGroupFile(""), sMatrixName(""), fillType(0),
                      sErrFile( "errors.txt" )
                      { algo[0] = ""; algo[1] = ""; }
}   ProgOptions;

typedef struct image_options {
	//string  image_pref;            
	int     iW       ;             // width.  Default 640 
	int     iH       ;             // height. Default 480 
	int     dX1      ;             // left   space; used by internal funcs 
	int     dX2      ;             // right  space 
	int     dY1      ;             // top    space 
	int     dY2      ;             // bottom space 
	int     gapColor ;             // color index
	int     unpairColor ;          // color index;  'single/mismatch'
    int     type     ;    // **reserved // 0 -std, 1 - max in col, 2 - selected gr, 3 - most gr, 4 - least gr
                          // 11 - stat
    int     offType  ;    // offset scale type; 0 - from 0 (abs), 1 - from min (relative)
    
    int     WLwidth;
    int     WLtype;       // 0 - bottom, 1 - inside; 2 - top (reseeved)

    // -----------------------------------------------------------------------------
    // some really unused - to satisfy AligraImage(), another for type = 22.
    int     fillType ;    // 0 - full (unweighted), 6 - weighted middle, 7 - very wei middle
    double              WeightLineMax;         //   for weighted filling
    double              Weight;                //
    double              WeightLineScaleMax;    //   for col values (incl. comparison picture)
                                               //   i.e. max of two aligns max for comp.pic.
                                               //   and otherwise eq. WeightLineMax. 
                                               //   Set inside AligraImage().

    vColGroupDef          *pColGroup;
    vColGroupDef          *pColGroup2;         // type = 22
    vGroupValueDef        *pGroupValue;
    vSelGroupDef          *pSelGroup;
    int     splitBlockSize;  // -1 - n/a, 0 - auto, 50, 60, 80, 100, 150 ... 500
    // end of AligraImage() variables
    // ---------------------------------------------------------------------------

    image_options( ) : iW(640), iH(480), type(0), fillType(0)
                       {}
    //image_options( ProgOptions& cfg) : iW(cfg.iW), iH(cfg.iH), type(cfg.vType[0]),
    //               pSelGroup( &cfg.vSel ), WLwidth(10), WLtype(0)     
    image_options( ProgOptions& cfg) : iW(cfg.iW), iH(cfg.iH), 
                   WLwidth(10), WLtype(0), offType(cfg.offType),
                   fillType(cfg.fillType), splitBlockSize(-1)
                       {}

}   ImageOptions;

typedef struct alignment {
    vSNamesDef          *pSNames;
    vOffsDef            *pOffs;
    vAliDef             *pAli;            // for f_outFragmentWL

    //+string              *pWeightLine;  // input WL *NIK rejected
    //+string              *pWLlegends;   // input WL *NIK rejected
    double              *WeightLineSums;  // [] ; computed w matrix
    double              WeightLineMax;    //
    vIntDef             *pMisPos;         // Mismatch (one) Positions (if any)

    // -----------------------------------------------------------------
    // Block to satisfy AligraImage() for type = 22
    vAliDef             *pAli2;           // Aligns w gaps   - for comp pic (type=22)
    double              *WeightLineSums2;  // for comparison picture (type = 22)
    double              WeightLineMax2;    //
    int                 WeightLineLen2;    // length maybe not eq to the 1st WL;
                                           // Maybe set less/greater/equal.

    vGroupDef           *pGroup;          // Group letters
    vGroupValueDef      *pGroupValue;     // reserved (computed and not used so far)
    vColGroupDef        *pColGroup;       // Dominant groups at columns
    // ------------------------------------------------------------------

    alignment() // : pWeightLine(0), pWLlegends(0) 
                  {}
} Align;

#include <alirow_iter.h>
#include <alico_image.h>

#include <aligets.h>

// ------------------------------------------------------------------
/*****
int ChoosParam( char par, char *LetPar, int *IntPar )   {
int num = strlen( LetPar );

        char Par = toupper( par );
        for( int i = 0; i < num; i++ )  {
                if( Par == LetPar[ i ] )
                        return IntPar[ i ];
        }
        // ******** THERE IS DEFAULT VALUE ***********
        return  IntPar[ 0 ];
}
*****/

int     ReadProgOptions( ProgOptions& cfg, char *file_name )      {
	    
	    char  *cEndPos;

        char buf[ _OPTIONS_STR_LEN ];
        FILE *f = fopen( file_name, "rt" );
        if( !f )        return 1;

        while( fgets( buf, _OPTIONS_STR_LEN, f ) )       {
                // from 1 to the first != ' '
			    //+cout << buf ;
				if( strlen( buf ) < 2 ) continue;   // <----------------

                char *par = buf + 2 + strspn( buf + 2, " " );
                switch( buf[1] )        {  // First symbol is "-"

               /***
                case   'h':
                        cfg.homology_only = ChoosParam( *par, LetParHom, IntParHom );
                        break;
					 case   'n':
								cfg.Naway = ( toupper( *par ) == 'Y' ) ? 1 : 0;
                        break;
					 case   'l':
								cfg.log_scale = ChoosParam( *par, LetParSca, IntParSca );
                        break;
                case   'a':
                        cfg.algorithm = ChoosParam( *par, LetParAlg, IntParAlg );
                        break;
                case   'f':
                        cfg.from = atoi( par );
                        break;
                case   't':
                        cfg.to = atoi( par );
                        break;
                case   'g':
                        cfg.gap_penalty = atoi( par );
                        break;

                // Tree Format to Output
                case   'o':
						if( strncmp( par, "PHYLIP", 5 ) == 0 ) {
                            tf.PHYLIP = 1;
                        } else if( strncmp( par, "PHYLMULTI", 5 ) == 0 ) {
                            tf.PHYLMULTI = 1;
                        } else if( strncmp( par, "NEXUS", 5 ) == 0 ) {
                            tf.NEXUS = 1;
                        } else if( strncmp( par, "NEXTREE", 5 ) == 0 ) {
                            tf.NEXTREE = 1;
                        } else if( strncmp( par, "HENNIG", 5 ) == 0 ) {
                            tf.HENNIG = 1;
                        } 
                        cout << " - Tree  Format: " << par; // << endl;
                        break;

                case   'i':
					    cEndPos = strpbrk( par, "\n\r" );
						if( cEndPos != 0 )
							*cEndPos = '\0';

								cfg.image_pref = par;           // string
						cout << " - Image   Pref: " << cfg.image_pref << endl;
								break;
				// ************** Changed by VNik ************* 02.05.2001
				// -u - removed (?)
				case   'C':
						if( strncmp( par, "Disable", 7 ) == 0 || strncmp( par, "DISABLE", 7 ) == 0  ) {
							cfg.unroot_disable = 1;
							cout << " - " << par+1 << endl;
						}
						cEndPos = strpbrk( par, "\n\r" );
						if( cEndPos != 0 )
							*cEndPos = '\0';

						sMM_Factor[0] = par; // string
						cout << " - Cluster Maxmin Factor: " << sMM_Factor[0] << endl;
						break;

				case   'T':
						if( strncmp( par, "Disable", 7 ) == 0 || strncmp( par, "DISABLE", 7 ) == 0  ) {
							cfg.unroot_disable = 1;
							cout << " - " << par+1 << endl;
						}
						cEndPos = strpbrk( par, "\n\r" );
						if( cEndPos != 0 )
							*cEndPos = '\0';

						sMM_Factor[1] = par;
						cout << " - Topo Maxmin Factor:    " << sMM_Factor[1] << endl;
						break;

				// ******************* Restoring -u :
                case   'u':
                        if( strncmp( par, "Disable", 7 ) == 0 || strncmp( par, "DISABLE", 7 ) == 0  ) {
                            //tf.UNROOT_DISABLE = 1;
							cfg.unroot_disable = 1;
						    cout << " - " << par+1 << endl;
							break;                          // *** inserted
						}
					    cEndPos = strpbrk( par, "\n\r" );
						if( cEndPos != 0 )
							*cEndPos = '\0';

                        if( strncmp( par, "Cluster Maxmin Factor ", 22 ) == 0 ) {
							sMM_Factor[0] = par + 21 + strspn( par + 21, " " ); // string
						    //cout << " - " << par << endl;
						    cout << " - Cluster Maxmin Factor " << sMM_Factor[0] << endl;
						}
                        if( strncmp( par, "Topo Maxmin Factor ", 19 ) == 0 ) {
							sMM_Factor[1] = par + 18 + strspn( par + 18, " " );
						    //cout << " - " << par << endl;
						    cout << " - Topo Maxmin Factor " << sMM_Factor[1] << endl;
						}

						break;
                // ****************************** end of -u
            ***/

				case   'W':
						cfg.iW = atoi( par );
						if( cfg.iW < _IMAGE_MIN_WIDTH ) {   // Check limits
							cfg.iW = _IMAGE_MIN_WIDTH;
						}
						if( cfg.iW > _IMAGE_MAX_WIDTH ) {
							cfg.iW = _IMAGE_MAX_WIDTH;
						}
						cout << " - Image  Width: " << cfg.iW << endl;
								break;
				case   'H':
						cfg.iH = atoi( par );
						//if( cfg.iH < _IMAGE_MIN_HEIGHT ) {   // Check limits
						//	cfg.iH = _IMAGE_MIN_HEIGHT;
						//}
						if( cfg.iH > _IMAGE_MAX_HEIGHT ) {
							cfg.iH = _IMAGE_MAX_HEIGHT;
						}
						cout << " - Image Height: " << cfg.iH << endl;
								break;
				case   'M':
						cfg.offMisMode = atoi( par );
						cout << " - Mismtch Mode: " << cfg.offMisMode << endl;
								break;
				case   'o':
						cfg.offType = atoi( par );
						cout << " -  Offset Type: " << cfg.offType << endl;
								break;
                case   'i':
					    cEndPos = strpbrk( par, "\n\r" );
						if( cEndPos != 0 )
							*cEndPos = '\0';

								cfg.image_pref = par;           // string
						cout << " - Image   Pref: " << cfg.image_pref << endl;
								break;
                case   'm':
					    cEndPos = strpbrk( par, "\n\r" );
						if( cEndPos != 0 )
							*cEndPos = '\0';

								cfg.sMatrixName = par;           // string
						cout << " - Matrix  File: " << cfg.sMatrixName << endl;
								break;
                case   'e':
					    cEndPos = strpbrk( par, "\n\r" );
						if( cEndPos != 0 )
							*cEndPos = '\0';

								cfg.sMismatchFile = par;           // string
						cout << " - Mismat. File: " << cfg.sMismatchFile << endl;
								break;
                case   'E':
					    cEndPos = strpbrk( par, "\n\r" );
						if( cEndPos != 0 )
							*cEndPos = '\0';

								cfg.sErrFile = par;           // string
						cout << " - Error   File: " << cfg.sErrFile << endl;
								break;
                case   'g':                                     // for type = 22  
					    cEndPos = strpbrk( par, "\n\r" );
						if( cEndPos != 0 )
							*cEndPos = '\0';

								cfg.sGroupFile = par;           // string
						cout << " - Group   File: " << cfg.sGroupFile << endl;
								break;
				case   't':                                  // type
						//cfg.type = atoi( par );
                        int ty;
						ty = atoi( par );
                        (cfg.vType).push_back(ty);
						cout << " - Image   Type: " << ty << endl;
								break;
				case   'f':                                  // type
                        
						cfg.fillType  = atoi( par );
                        
						cout << " - Filling Type: " << cfg.fillType << endl;
								break;
            /**
				case   's':                                  // Selected groups
						//cfg.type = atoi( par );
                        int iSel;
						iSel = atoi( par );
                        (cfg.vSel).push_back(iSel);
						cout << " - Sel.   Group: " << iSel << endl;
            **/
								break;
				case   'T':                                  // title
					    cEndPos = strpbrk( par, "\n\r" );
						if( cEndPos != 0 )
							*cEndPos = '\0';
						cfg.title = par;
						cout << " - Image  Title: " << cfg.title << endl;
								break;
				case   'A':                                  // algo1 name
					    cEndPos = strpbrk( par, "\n\r" );
						if( cEndPos != 0 )
							*cEndPos = '\0';
						cfg.algo[0] = par;
						cout << " -1st Algorithm: " << cfg.algo[0] << endl;
								break;
				case   'B':                                  // algo2 name
					    cEndPos = strpbrk( par, "\n\r" );
						if( cEndPos != 0 )
							*cEndPos = '\0';
						cfg.algo[1] = par;
						cout << " -2nd Algorithm: " << cfg.algo[1] << endl;
								break;
				}
		  }
		  cout << "----------------------------------------------" << endl << endl;

          //if( (cfg.vSel).size() > 0 )
          //    (cfg.vType).push_back( 2 );          // add type "sel groups"
          return 0;
}

/*****
int   f_getAliOffs_(  char *fAli, vSNamesDef& vSNames1, vAliDef& vAli1,
                     vOffsDef& vOffs1, string& sWL ) {

    // Output:
    //      vOffs1: gaps and 'e' only
    //       vAli1: w/o gaps 

    mOffDef                                     mOff;
    mOffDef::iterator                           iiO;
    int                 i, iP;

    //FILE *afile = fopen( argv[ 1 ], "rt" );
	FILE *afile = fopen( fAli, "rt" );
	if( !afile )    {
					 cout << "\n Can't open align file " << fAli << endl;
					 return 2;
	}
    f_getAlign( afile, vSNames1, vAli1, sWL );
    getch();


    // Strip alignments off gaps & create offset map w only gaps.
    int         iSize; 
    for( i = 0; i < vAli1.size(); ++i ) {
        iP = (vAli1[i]).rfind( '-' );
        cout << "*** " << vAli1[i] << ", iP- = " << iP << endl;

        mOff.clear();
        //j = 0;
        iSize = (vAli1[i]).length();
        while( (iP = (vAli1[i]).rfind( '-' ) ) > -1 ) {
        // - ***err! while( (iP = (vAli1[i]).find( '-' ) ) > -1 ) {  //* next pos changed ! * err
            mOff[iP]  = OffsetDescr( 'g', 1 );
            (vAli1[i]).erase( iP, 1 );
        }
        mOff[iSize]  = OffsetDescr( 'e' );

        vOffs1.push_back( mOff );
    }
    getch();


    // *tst - display alignments w/o gaps and offsets including 'e'.
    for( i = 0; i < vAli1.size(); ++i ) {
        cout <<  endl << "--- " << vAli1[i] << endl;
        cout <<          "  - Offs: " ;
        for( iiO = (vOffs1[i]).begin(); iiO != (vOffs1[i]).end(); ++iiO ) {
            cout << (*iiO).first << "  " ;
        }
        cout << endl;
    }
    getch();

    fclose( afile );

    return 0;
}
*****/

//int   f_getAliOffs(  char *fAli, vSNamesDef& vSNames1, vAliDef& vSeq1, vAliDef& vAli1,
int   f_getAliOffs(  FILE *afile, vSNamesDef& vSNames1, vAliDef& vSeq1, vAliDef& vAli1,
                     vOffsDef& vOffs1, string& sWL ) {

    // Output:
    //      vOffs1: gaps and 'e' only
    //       vAli1: w/o gaps 

    mOffDef                                     mOff;
    mOffDef::iterator                           iiO;
    int                 i, iP;

	//FILE *afile = fopen( fAli, "rt" );
	//if( !afile )    {
	//				 cout << "\n Can't open align file " << fAli << endl;
	//				 return 2;
	//}
    f_getAlign( afile, vSNames1, vAli1, sWL );
    _GETCH;

    vSeq1 = vAli1;         // **************

    // Strip alignments off gaps & (-- create offset map w only gaps).
    int         iSize; 
    for( i = 0; i < vSeq1.size(); ++i ) {
        iP = (vSeq1[i]).rfind( '-' );
#ifdef _TRACE
        cout << "*** " << vSeq1[i] << ", iP- = " << iP << endl;
#endif
        mOff.clear();
        //j = 0;
        iSize = (vSeq1[i]).length();
        while( (iP = (vSeq1[i]).rfind( '-' ) ) > -1 ) {
        // - ***err! while( (iP = (vAli1[i]).find( '-' ) ) > -1 ) {  //* next pos changed ! * err
            //++ mOff[iP]  = OffsetDescr( 'g', 1 );
            (vSeq1[i]).erase( iP, 1 );
        }
        mOff[iSize]  = OffsetDescr( 'e' );

        vOffs1.push_back( mOff );
    }
    _GETCH;


    // *tst - display alignments w/o gaps and offsets including 'e'. Plus w gaps.
    for( i = 0; i < vSeq1.size(); ++i ) {
#ifdef _TRACE
        cout <<  endl << "--- " << vSeq1[i] << endl;
        //cout <<  endl << " -a " << vAli1[i] << endl;
#endif
      /**
        cout <<          "  - Offs: " ;
        for( iiO = (vOffs1[i]).begin(); iiO != (vOffs1[i]).end(); ++iiO ) {
            cout << (*iiO).first << "  " ;
        }
        cout << endl;
      **/
    }
    _GETCH;

    fclose( afile );

    return 0;
}


/*****

// 'g' of self, and of contrepart as 'd', and 'e'.
//  And 'u' (if any, and w/o 'd').
void    f_setOneSeqDelta( mOffDef&   mOff ) {

    int             i, iSize=0, iSum=0, iP, iP0;
    char            ty, ty0;

    mOffDef::iterator                           iiO, iiO0;
    OffsetDescr                                 oDescr, oD0, oD;

    bool                                        bUnpaired = false;

    iiO     = mOff.begin();               // at least 'e'
    oD0     = (*iiO).second;
    ty0     = oD0.ty;
    if( ty0 == 'u' || ty0 == 'm' ) {            
       bUnpaired = true;
       iiO++;
       oD0     = (*iiO).second;
       ty0     = oD0.ty;
    }

    iP0     = (*iiO).first;
    iSum    = oD0.ddx;
   
    // 'g' -> 0
    if( ty0 == 'g' ) {            // 'd' -1 : left unchanged !
        oD0.ddx = 0;              // +1 already added to iSum
        (*iiO).second = oD0;
    }

    cout << "*** Set One Seq Delta" << endl;

    iiO0 = iiO;
    for( ++iiO; iiO != mOff.end(); ++iiO ) {

        oD     = (*iiO).second;
        iP     = (*iiO).first;
        ty     = oD.ty;

        // to be done some other day...

        //if(  ty0 == 'm' ) {            
        //   bUnpaired = true;
        //   continue;
        //}


        // Break of g/d
        if( iP > (iP0 + 1 ) ) {   
        //} else {

            // Insert after gap
            if( ty0 == 'g' ) {
                oD0.ddx  = iSum;
                oD0.ty   = 'd';
                if( iSum == 0 )
                    oD0.ty   = '0';
                if( bUnpaired ) {             // ??????????
                    oD0.ty   = '0';
                    oD0.ddx  = 0;
                }

                mOff[iP0+1] = oD0;    // insert before iP ***** iiO is valid 
                cout << " -i- ty/pos/ddx = " << oD0.ty << ":" << iP0 + 1 << "/" << iSum << endl;

            } else if( ty0 == 'd' ) {
                oD0.ddx  = iSum;
                oD0.ty   = 'd';
                if( iSum == 0 )
                    oD0.ty   = '0';
                if( iSum < 0 ) {
                    mOff[iP0] = oD0;    // change (if any) ddx of same pos  
                    cout << " -ch/d/P0- ty/pos/ddx = " << oD0.ty << ":" << iP0 << "/" << iSum << endl;
                } else {
                    //mOff.erase( iiO0 );   // erase pos iP0
                    //cout << " -e  pos << iP0 << endl;

                    mOff[iP0+1] = oD0;    // change (if any) ddx of same pos  
                    cout << " -i/d/P0+1- ty/pos/ddx = " << oD0.ty << ":" << iP0+1 << "/" << iSum << endl;
                }

            }
        }

        //// g/d without space
        //if( iP == (iP0 + 1 ) ) {    // gggd      NOT gggddgdd
            iSum    += oD.ddx;
            if( ty == 'g' ) {
                oD.ddx = 0;
                (*iiO).second = oD;
                cout << " -g- ty/pos/ddx = " << oD.ty << ":" << iP << "/" << 0 << endl;
            } else if( ty == 'd' ) {
                //oD.ddx = iSum;
                oD.ddx = iSum + 1;
                if( (iSum + 1) == 0 ) 
                    oD.ty = '0';
                (*iiO).second = oD;
                cout << " -d/0+ ty/pos/ddx = " << oD.ty << ":" << iP << "/" << iSum + 1 << endl;

            }
        //}

        ty0  = ty;
        iP0  = iP;
        iiO0 = iiO;
        //getch();

    }

    // Remove next doubles
    iiO  = mOff.begin();
    iiO0 = iiO;
    iP0  = (*iiO0).first;
    oD0  = (*iiO0).second;

    for( ++iiO; iiO != mOff.end();  ) {
        oD  =  (*iiO).second;
        if( oD0.ty == oD.ty ) {
            //if( oD.ty == 'g' && (iP0 + 1) >= (*iiO).first ) {
            if( oD.ty == 'g' ) {
                mOff.erase( iiO++ );  // erased iiO

            } else if( oD.ty == 'd' && (oD0.ddx == oD.ddx) ) {
                mOff.erase( iiO++ );  // erased iiO
            } else {
                //iiO0 = iiO;
                iP0 = (*iiO).first;
                ++iiO;
            }
        } else {
            //iiO0 = iiO;
            iP0 = (*iiO).first;
            ++iiO;
        }
    }

}

void    f_calcDelta_( vOffsDef& vOffs1, vOffsDef& vOffs2 , vRefDef& vRef12 )

//                     mSNamesDef&  mSN1 ) // 1st arg changed, 2nd not changed
{
    int             i, iSize=0;
    //mOffDef                                     mOff;
    mOffDef::iterator                           iiO;
    OffsetDescr                                 oDescr, oD1;


    // Add gaps of the contrepart
    for( i = 0; i < vOffs1.size(); ++i ) {

        mOffDef&            mOff = vOffs1[i];
        int                 iP   = vRef12[i];

        // Unpaired - no contrepart
        if(  iP < 0 ) {
            if( iP == -2 ) {
                mOff[-1]  = OffsetDescr( 'm' );
                cout << "*** Mismatch " << i << endl;
                getch();
            } else {
                mOff[-1]  = OffsetDescr( 'u' );
                cout << "*** Unpaired " << i << endl;
            }

            f_setOneSeqDelta( mOff );
            continue;                        // <------------
        }

        // Paired
        iiO = mOff.end();
        iiO--;
        iSize = (*iiO).first - 1;    // The last item must be 'e'
#ifdef _TRACE
        cout << "--- Pair " << i << ", " << iP << ": " ; //endl;
#endif
        // Insert gaps of the contrepart
        //for( iiO = (vOffs2[iP]).begin(); iiO != (vOffs2[iP]).end(); ++iiO ) {
        for( iiO = (vOffs2[iP]).begin(); iiO != (vOffs2[iP]).end();  ) {
            int       iG00, iG0, iG = (*iiO).first;
            oDescr    = (*iiO).second;

            if( iG > iSize )           break;     // -----------> contrepart longer ??

            if( oDescr.ty != 'g' )   {
                ++iiO;
                continue;  // <------------
            }

            // Get target:
            oD1       = mOff[iG];   // ddx = +1 if found, 0(implicit insert) otherwise
            if( oD1.ty != 'g' )
                oD1.ty  = 'd';      // if not found in target mOff

            --oD1.ddx;              // 0 or -1

            if( oD1.ty == 'g' ) {
                mOff[iG]  = oD1;        // replace with ddx = 0
#ifdef _TRACE
                cout << iG << "/" << oD1.ddx; 
#endif
                ++iiO;

            // Get all nearby contrepart gaps until a target gap if any:
            } else {
                iG00 = iG0 = iG;
                ++iiO;
                while( iiO != (vOffs2[iP]).end() ) {
                    iG      = (*iiO).first;
                    if( iG > iSize )
                        break;       // ----------> contrepart longer

                    oDescr  = (*iiO).second;
                    if( iG > (iG0 + 1 ) ) break;
                    if( oDescr.ty != 'g' ) break;

                    --oD1.ddx;
                    iG0 = iG;
                    ++iiO;
                }

                // Insert 'd' (really replace implicit defaults):
                mOff[iG00]  = oD1;      
#ifdef _TRACE
                cout << iG00 << "/" << oD1.ddx; 
#endif                
            }
        }
        cout << endl;

        // Set offsets
        f_setOneSeqDelta( mOff );

    }

    getch();
}
*****/

void    f_calcDelta( vOffsDef& vOffs1, vOffsDef& vOffs2 , vRefDef& vRef12,
                     vAliDef&  vAli1, vAliDef&  vAli2 )

//                     mSNamesDef&  mSN1 ) // 1st arg changed, 2nd not changed
{
    int             i, iSize=0; // *tst , j ,iWk, iWk2;
    //mOffDef                                     mOff;
    mOffDef::iterator                           iiO;
    OffsetDescr                                 oDescr, oD1, oD2;


    // Add offs (except gaps) //of the contrepart
    for( i = 0; i < vOffs1.size(); ++i ) {

        mOffDef&            mOff = vOffs1[i];             // +- ?  *** 28.06.2001
        //++mOffDef             mOff = vOffs1[i];
        int                 iP   = vRef12[i];
        //bool                bUnpair = false;


        AlignRowIterator   iiA( &( vAli1[i] ) ), iiB(  &( vAli2[iP] ) );
        int                iOff0, iOff;
        int                iLen1=(vAli1[i]).length(), iLen2 = (vAli2[iP]).length();

        // Unpaired - no contrepart
        if(  iP < 0 ) {
            if( iP == -2 ) {
                mOff[-1]  = OffsetDescr( 'm' );
                cout << "*** Mismatch " << i << endl;
                _GETCH;
            } else {
                mOff[-1]  = OffsetDescr( 'u' );
                cout << "*** Unpaired " << i << endl;
            }

            //bUnpair = true;
            for( iiA.begin(); iiA.align_pos() < iLen1; ++iiA ) {
                if( iiA.gaps_before() > 0 ) {
                    mOff[*iiA - iiA.gaps_before()] = OffsetDescr( 'g' );
                    mOff[*iiA] = OffsetDescr( '0' );
                }
            }
            // Gaps at end
            if( iiA.gaps_before() > 0 ) {
                mOff[*iiA - iiA.gaps_before()] = OffsetDescr( 'g' );
            }
            continue;                        // <------------
        }

        // Paired
        iiO = mOff.end();
        iiO--;
        iSize = (*iiO).first - 1;    // The last item must be 'e'
#ifdef _TRACE
        cout << "--- Pair " << i << ", " << iP << ": " ; //endl;
#endif

        // insert offsets  
        //iiA0  = iiA.begin();
        //iiB0  = iiB.begin();
        iiA.begin();
        iiB.begin();
        iOff0 = -iLen1 - 1;   // to be not eq the next offset

        // *tst for( j=0; (iiA.align_pos() < iLen1) && (iiB.align_pos() < iLen2); ++iiA, ++iiB, ++j ) {
        for( ; (iiA.align_pos() < iLen1) && (iiB.align_pos() < iLen2); ++iiA, ++iiB ) {
            iOff  = *iiA - *iiB;

            //*tst
            //++iWk  = *iiA;   // *tst iWk, iWk2, oD2

            if( iOff != iOff0 ) {

              if( iiA.gaps_before() > 0 ) {
                    // *tst
                    //++iWk2 = iiA.gaps_before();
                    //cout << "*** GapsBefore: *iiA/gaps =" << iWk << "/" << iWk2 << endl;
                    //getch();
                    //++oD2 = mOff[*iiA - iiA.gaps_before()];
                    mOff[*iiA - iiA.gaps_before()] = OffsetDescr( 'g' );
                    //mOff.insert( mOffDef::value_type(*iiA - iiA.gaps_before(), OffsetDescr( 'g' )) );
              }

              if( iOff == 0 ) {
                mOff[*iiA] = OffsetDescr( '0' );

#ifdef _TRACE
                cout << "  " << *iiA << ":0" ;
#endif
              } else {
                mOff[*iiA] = OffsetDescr( 'd', iOff );
#ifdef _TRACE
                cout << "  " << *iiA << ":" << iOff ;
#endif
              }
            }
            iOff0 = iOff;
        }
        // Gaps at end:
        if( iiA.gaps_before() > 0 ) {
                    mOff[*iiA - iiA.gaps_before()] = OffsetDescr( 'g' );
        }

#ifdef _TRACE
        cout << endl;
#endif
        // Save offsets:
        //++vOffs1[i] = mOff;                // *** 28.06.2001
    }

#ifdef _TRACE
    _GETCH;
#endif
}

long  f_strCount( string& s, char c ) {
    long   i,  iL = s.length(), cou = 0;

    for( i = 0; i < iL; ++i ) {
       if( s[i] == c ) ++cou;
    }
    return cou;
}

/***
// confirm input WL                                  *** NIK rejected - 15.06.2001
void    f_confWL( string& sWL1, string& sWLleg1 ) {
#define     _SW_LEGEND_GENEBEE     " .+*"
#define     _SW_LEGEND_CLUSTAL     " .:*"
    long    iCou[5], iTotal = 0;  // 0 - blanl, 1 - point, 2 - :, 3 - + , 4 - *
    
    iCou[0] = f_strCount( sWL1, ' ' ); iTotal  = iCou[0];
    iCou[1] = f_strCount( sWL1, '.' ); iTotal += iCou[1];
    iCou[2] = f_strCount( sWL1, ':' ); iTotal += iCou[2];
    iCou[3] = f_strCount( sWL1, '+' ); iTotal += iCou[3];
    iCou[4] = f_strCount( sWL1, '*' ); iTotal += iCou[4];
    cout << " -- Count: iTotal = " << iTotal << "(" << iTotal - iCou[0] 
         << "), len = " << sWL1.length() << endl;

    // OK: not only blanks and enough others
    if( iCou[0] < iTotal && iTotal * 2 > sWL1.length() ) {
        if( iCou[3] > 0 ) 
            sWLleg1 = _SW_LEGEND_GENEBEE;
        else if( iCou[2] > 0 )
            sWLleg1 = _SW_LEGEND_CLUSTAL;
        else
            sWLleg1 = _SW_LEGEND_GENEBEE;

        cout << "  - WL confirmed" << endl;

    // NOT confirmed:
    } else {
        sWL1 = "";
        sWLleg1 = "";
        cout << "  * WL NOT confirmed" << endl;
    }
}
***/

int     f_posAlign( string& a, int p0 ) {   // p0 - pos of seq, 
    int     i, j = 0;  // j - align pos to be found
    for( i = 0; i < a.length(); ++i ) {
        if( a[i] == '-' ) continue;   // <------------- gap
        ++j;
        if( j >= p0 )
            break;       // ------> found
    }
    return i + 1;
}

// s1, s2 - seqs (w/o gaps), a1, a2 - align seqs (with gaps)
int    f_outMismatch( FILE *f, int iS1, int iS2, string& sName, string& s1, string& s2,
                       string& a1, string& a2, int& iMP1, int& iMP2 ) {
    int         i, iP = -1, iL = s1.length();
    if( s2.length() < iL )
        iL = s2.length();

    // Find 1st letters:
    for( i = 0; i < iL; ++i ) {
        if( s1[i] != s2[i] ) {
            iP = i;
            break;
        }
    }

    if( iP > -1 ) {
        iMP1 = f_posAlign( a1, iP );
        iMP2 = f_posAlign( a2, iP );
        fprintf( f, "%s (%d,%d) - position in sequence: %d,  letter <font color=red>%c</font> (pos %d) versus <font color=red>%c</font> (pos %d).\n",
             sName.c_str(), iS1, iS2, iP, s1[iP], iMP1, s2[iP], iMP2 ); 


    // No different leters
    } else {
        iMP1 = iMP2 = -1;
        if( s1.length() == s2.length() ) {                          // *** ?
            fprintf( f, "%s (%d,%d) - mismatch of the two sequences.\n",
                 sName.c_str(), iS1, iS2 );
            iP   = -2;  // *undef
        } else {
            fprintf( f, "%s (%d,%d) - different lenghts: <font color=red>%d</font> versus <font color=red>%d</font>.\n",
                 sName.c_str(), iS1, iS2, s1.length(), s2.length() );
            iP = -1;  // dif. len
            if( s1.length() > s2.length() )
                iMP1 = s2.length() + 1;
            else if( s1.length() < s2.length() )
                iMP2 = s1.length() + 1;
        }
    }
    return iP;
}


int main( int argc, char** argv )   {

	if( argc < 4 )	{
		printf( "Compare Two Alignments\n " );
		printf( "USAGE: alicomp.exe optName alignName1 alignName2 \n" );
		printf( "      optName  - fileName of alignment options\n" );
		printf( "    alignName  - fileName of input alignment\n" );
		return	1;
	}

	// *****************  Update parameters  *****************
	ProgOptions          cfg; 
	if( ReadProgOptions( cfg, argv[ 1 ] ) )   {
			 cout << "\n Can't read params file: " << argv[ 1 ]  << endl;
             return 1;
	}
    
    ImageOptions   image_cfg( cfg );

    //vector< string, alloc >                                  vSNames;
    //vector< map< int, off_de, less< int >, alloc >, alloc >  vOffs;
    //map< int, off_de, less< int >, alloc >                   mOff;
    vAliDef                                     vAli1, vAli2, vSeq1, vSeq2;
    vSNamesDef                                  vSNames, vSNames1, vSNames2;
    vOffsDef                                    vOffs, vOffs1,  vOffs01, vOffs2;
    mOffDef                                     mOff;
    mOffDef::iterator                           iiO;
    OffsetDescr                                 oDescr; //, oD1;
    mSNamesDef                                  mSN, mSN1, mSN2;
    
    int                 i, j, iP;

    // tst ***************************************************************
    char   *sTit   = "L-lactate dehydrogenase sequences";
    char   *sTit5n = "Sample from Tree Help";
    char   *sAlgo1 = "GeneBee";
    char   *sAlgo2 = "ClustalW";
    char   *sAlgoAE2 = "AEGOR 2";
    char   *sAlgoAE3 = "AEGOR 3";

    Align             oAli, oAli1, oAli2;
    string            sWL, sWLleg, sWL1, sWLleg1, sWL2, sWLleg2;

// **
	FILE *af1 = fopen( argv[ 2 ], "rt" );
	if( !af1 )    {
					 cout << "\n Can't open align file " << argv[ 2 ] << endl;
					 return 2;
	}
	FILE *af2 = fopen( argv[ 3 ], "rt" );
	if( !af2 )    {
					 cout << "\n Can't open align file " << argv[ 3 ] << endl;
					 return 2;
	}
    //f_getAliOffs(  argv[ 2 ], vSNames1, vSeq1, vAli1, vOffs1, sWL1 );
    //f_getAliOffs(  argv[ 3 ], vSNames2, vSeq2, vAli2, vOffs2, sWL2 );
    f_getAliOffs(  af1, vSNames1, vSeq1, vAli1, vOffs1, sWL1 );
    f_getAliOffs(  af2, vSNames2, vSeq2, vAli2, vOffs2, sWL2 );
    // vOffs's contain just gaps w 'ddx = +1' and end-mark 'e'

     // Confirm WL's & type of symbol set
    
    //+f_confWL( sWL1, sWLleg1 );  // *NIK rejected
    //+f_confWL( sWL2, sWLleg2 );
    //+getch();


    int          nSeq1 = vAli1.size(), iSeqLen1 = (vAli1[0]).length();
    double       *sum1 = new double[ iSeqLen1 ];
    int          nSeq2 = vAli2.size(), iSeqLen2 = (vAli2[0]).length();
    double       *sum2 = new double[ iSeqLen2 ];

    int          molType = f_getMolType( vAli1 );

    // matrix file - AUTO
    if( cfg.sMatrixName == "" ) {
        cfg.sMatrixName = "blosum62.matrix";
        if( molType == 1 )
            cfg.sMatrixName = "dnarna.matrix";
    }

    // Check out molType & matrix consistency:             *** 22.08.2001
    FILE                    *fErr;
    string                  sMT_ERR = "";
    if( molType == 1 && cfg.sMatrixName != string("dnarna.matrix") ) {
      sMT_ERR = "The sequences are found to be nucleotide. The matrix \"" + 
                string( cfg.sMatrixName ) + "\" is protein.\nResult may be inconsistent.";

    } else if ( molType != 1 && cfg.sMatrixName == "dnarna.matrix" ) {
      sMT_ERR = "The sequences are found to be protein. The matrix \"" + 
                string( cfg.sMatrixName ) + "\" is nucleotide.\nResult may be inconsistent.";

    }

    if( !sMT_ERR.empty() ) {
        fErr = fopen( (cfg.sErrFile).c_str(), "wt" );
        fprintf( fErr, "%s\n", sMT_ERR.c_str() );
        fclose( fErr );  // ?
    }


    //oAli1.WeightLineMax  =  f_setWeights( vAli1, "blosum62.matrix", sum1 );
    //oAli2.WeightLineMax  =  f_setWeights( vAli2, "blosum62.matrix", sum2 );
    oAli1.WeightLineMax  =  f_setWeights( vAli1, (cfg.sMatrixName).c_str(), sum1 );
    oAli2.WeightLineMax  =  f_setWeights( vAli2, (cfg.sMatrixName).c_str(), sum2 );
    if( oAli2.WeightLineMax < 0 ) 
        return 3;   // ------------------------->                 // matrix ?
    
    //+getch();
    //_GETCH;
    //cout << " -- weight max = " << oAli1.WeightLineMax << ", " << oAli2.WeightLineMax << endl;

    vRefDef                 vRef12( vSNames1.size(), -1 );   // -1 - no ref
    vRefDef                 vRef21( vSNames2.size(), -1 );

    int                     iMismatch = 0, iWk;
    int                     iMP1, iMP2;            // mismatch pos
    FILE                    *fMismatch;

    // Mismatch Positions (one per seq at most)
    vIntDef                 vMisPos1( nSeq1, -1 );        // -1 - no mismatch 
    vIntDef                 vMisPos2( nSeq2, -1 );        // 

    // 1 <--> 2
    // mSN1/2 will contain pos of the contrepart or -1 otherwise:
    for( i = 0; i < vSNames2.size(); ++i ) {
        mSN.insert ( mSNamesDef::value_type( vSNames2[i], i + 1 ) );  // i + 1 *** !!!
        mSN2.insert( mSNamesDef::value_type( vSNames2[i], -1    ) );  // n/a
    }

    for( i = 0; i < vSNames1.size(); ++i ) {
        iP =  mSN[ vSNames1[i] ];    // if > 0  - okay
        mSN1.insert( mSNamesDef::value_type( vSNames1[i], iP - 1 ) );  // maybe -1 (n/a)

        // 1st found by Name in 2nd in pos iP:
        if( iP > 0 ) {
           
          // Match
          //+if( vAli1[i] == vAli2[iP - 1] ) {
          if( vSeq1[i] == vSeq2[iP - 1] ) {
            mSN2[ vSNames1[i] ] = i;     // set relationship if found
#ifdef _TRACE
            cout << " -- Rel: " << i << " <--> " << iP - 1 << endl;
#endif

            vRef12[i]      = iP - 1;
            vRef21[iP - 1] = i;

          // Mismatch
          } else {
            mSN1[ vSNames1[i] ] = -2;
            mSN2[ vSNames1[i] ] = -2;

            vRef12[i]      =  -2;
            vRef21[iP - 1] =  -2;

            cout << " ** Rel Mismatch: " << i << " <--> " << iP - 1 << endl;
            ++iMismatch;
            if( iMismatch == 1 ) {
                //fMismatch = fopen( "mismatch.txt", "wt" );
                fMismatch = fopen( (cfg.sMismatchFile).c_str(), "wt" );
            }
            //f_outMismatch( fMismatch, i, iP - 1, vSNames1[i], vAli1[i], vAli2[iP - 1] );
            iWk = f_outMismatch( fMismatch, i, iP - 1, vSNames1[i], vSeq1[i], vSeq2[iP - 1],
                           vAli1[i], vAli2[iP - 1], iMP1, iMP2 );

            if( iMP1 > -1 )
                vMisPos1[i] = iMP1;

            if( iMP2 > -1 )
                vMisPos2[iP - 1] = iMP2;



/*** ?                                                            // *** 28.06.2001
            // Only dif. lengths
            if( iWk == -1 ) {         // iWk >= 0 - dif. letters pos

                mSN2[ vSNames1[i] ] = i;     // set relationship if found
#ifdef _TRACE
                cout << " -- Rel (*dif. lengths): " << i << " <--> " << iP - 1 << endl;
#endif

                vRef12[i]      = iP - 1;
                vRef21[iP - 1] = i;
            }
***/
            // Mismatched still to set relationship
            if( cfg.offMisMode == 1 ) {
                mSN2[ vSNames1[i] ] = i;     // set relationship if found
#ifdef _TRACE
                cout << " -- Rel (*mismatched): " << i << " <--> " << iP - 1 << endl;
#endif
                vRef12[i]      = iP - 1;
                vRef21[iP - 1] = i;

            }

          }
        // Not found by Name - try compare vSeq's  //vAli's
        } else {
           cout << " ** Rel: " << i << " <--> *** Not Found by Name" << endl;
           iP = -1;
           for( j = 0; j < vAli2.size(); ++j ) {
               //+if( vAli1[i] == vAli2[j] ) {
               if( vSeq1[i] == vSeq2[j] ) {
                   iP = j;
                   break;
               }
           }
           if( iP > -1 ) {
               cout << "  * Rel: " << i << " <--> *** Still Found " << iP << endl;
               _GETCH;

               mSN2[ vSNames2[iP] ] = i;     // set relationship if found
               mSN1[ vSNames1[i ] ] = iP;     // set relationship if found

               vRef12[i]      = iP;
               vRef21[iP] = i;
           }

        }
    }
    if( iMismatch > 0 )
        fclose( fMismatch );


    // Add gaps of the contrepart  & then calc delta 
    vOffs01 = vOffs1;                        // Save
    f_calcDelta( vOffs1, vOffs2 , vRef12, vAli1, vAli2 ); //, mSN1 );  // 1st arg changed, 2nd not changed
    f_calcDelta( vOffs2, vOffs01, vRef21, vAli2, vAli1 ); //, mSN2 );  // 1st arg changed, 2nd not changed

    cout << endl << "*** 1st Offsets:" << endl;
    _GETCH;
    for( i = 0; i < vOffs1.size(); ++i ) {
        for( iiO = (vOffs1[i]).begin(); iiO != (vOffs1[i]).end(); ++iiO ) {
            oDescr = (*iiO).second;
            cout << "  " <<  (*iiO).first << ":'" << oDescr.ty << "':" << oDescr.ddx;
        }
        cout << endl;
    }
    cout << endl << "*** 2nd Offsets:" << endl;
    _GETCH;
    for( i = 0; i < vOffs2.size(); ++i ) {
        for( iiO = (vOffs2[i]).begin(); iiO != (vOffs2[i]).end(); ++iiO ) {
            oDescr = (*iiO).second;
            cout << "  " <<  (*iiO).first << ":'" << oDescr.ty << "':" << oDescr.ddx;
        }
        cout << endl;
    }
    //+getch();

   //oAli1.pWeightLine = &sWL1;
   //oAli1.pWLlegends  = &sWLleg1;

   oAli1.pSNames     = &vSNames1;
   oAli1.pOffs       = &vOffs1;

   //oAli2.pWeightLine = &sWL2;
   //oAli2.pWLlegends  = &sWLleg2;

   oAli2.pSNames     = &vSNames2;
   oAli2.pOffs       = &vOffs2;

   oAli1.WeightLineSums =  sum1;
   oAli2.WeightLineSums =  sum2;

   oAli1.pMisPos        =  &vMisPos1;
   oAli2.pMisPos        =  &vMisPos2;

   oAli1.pAli           =  &vAli1;         // for f_outFragmentWL() 
   oAli2.pAli           =  &vAli2;

    //+image_cfg.iH = 300; //250; // 150; 180; //300;   // 200  // 600
    //++AlicompImage( oAli1, image_cfg, "la10.png", sTit, sAlgo1, sAlgo2 );
    //++AlicompImage( oAli2, image_cfg, "la10b.png", sTit, sAlgo2, sAlgo1 );

    //+image_cfg.iH = 200;
    //AlicompImage( oAli1, image_cfg, "ftree3.png", sTit5n, sAlgo1, sAlgo2 );
    //AlicompImage( oAli2, image_cfg, "ftree3b.png", sTit5n, sAlgo2, sAlgo1 );

    //+image_cfg.iH = 200;

   image_cfg.WLwidth = 16;  //20;
   image_cfg.WLtype  = 0;  // 0 - bottom, 1 - inside

   if( cfg.algo[0] == "" ) 
       cfg.algo[0] =  "1st";
   if( cfg.algo[1] == "" ) 
       cfg.algo[1] =  "2nd";


   if( (cfg.vType).size() == 0 )
       (cfg.vType).push_back( 0 );   // defaults type  = 0

   // Auto height
   bool       bAutoHei = false; 
   int        iDefHei1 = 26 + 16 + 2  + nSeq1 * _IMAGE_DEFAULT_HEIGHT_DY +
                         8 + 16 + 16 + image_cfg.WLwidth + 12 + 12; 
   int        iDefHei2 = 26 + 16 + 2  + nSeq2 * _IMAGE_DEFAULT_HEIGHT_DY +
                         8 + 16 + 16 + image_cfg.WLwidth + 12 + 12; 
   if( iDefHei1 < _IMAGE_MIN_HEIGHT )
       iDefHei1 = _IMAGE_MIN_HEIGHT;
   if( iDefHei2 < _IMAGE_MIN_HEIGHT )
       iDefHei2 = _IMAGE_MIN_HEIGHT;
   if( image_cfg.iH <= 0 ) 
       bAutoHei = true;

   // Min Height DY
   int        iMinHei1 = 26 + 16 + 2  + nSeq1 * _IMAGE_MIN_HEIGHT_DY +
                         8 + 16 + 16 + image_cfg.WLwidth + 12 + 12; 
   int        iMinHei2 = 26 + 16 + 2  + nSeq2 * _IMAGE_MIN_HEIGHT_DY +
                         8 + 16 + 16 + image_cfg.WLwidth + 12 + 12; 
   if( iMinHei1 < _IMAGE_MIN_HEIGHT )
       iMinHei1 = _IMAGE_MIN_HEIGHT;
   if( iMinHei2 < _IMAGE_MIN_HEIGHT )
       iMinHei2 = _IMAGE_MIN_HEIGHT;

   // Save iH parm
   int     iH0 = image_cfg.iH, iH1, iH2;

   // Set iH1, iH2
   iH1 = iH2 = iH0;
   if( bAutoHei )
       iH1 = iDefHei1;
   else if( iH1 < iMinHei1 )
       iH1 = iMinHei1;

   if( bAutoHei )
       iH2 = iDefHei2;
   else if( iH2 < iMinHei2 )
       iH2 = iMinHei2;

   if( f_findType( cfg.vType, 0 ) ) {
        //cout << "*** Algo: " << cfg.algo[0] << "  " << cfg.algo[1] << endl;
        //+AlicompImage( oAli1, image_cfg, "aeg1.png", "", sAlgoAE2, sAlgoAE3 );
        //+AlicompImage( oAli2, image_cfg, "aeg1b.png", "", sAlgoAE3, sAlgoAE2 );
    
        image_cfg.iH = iH1;
        AlicompImage( oAli1, image_cfg, (cfg.image_pref+string("a.png")).c_str(), 
                 (cfg.title).c_str(), (cfg.algo[0]).c_str() , (cfg.algo[1]).c_str() );

        image_cfg.iH = iH2;
        AlicompImage( oAli2, image_cfg, (cfg.image_pref+string("b.png")).c_str(), 
                 (cfg.title).c_str(), (cfg.algo[1]).c_str() , (cfg.algo[0]).c_str() );
   }
   // ---------------------------------------------------------------------------------
   // type = 21/22 of AligraImage()
   //image_cfg.type = 21;  //22;      // see below

   if( !f_findType( cfg.vType, 21 ) && !f_findType( cfg.vType, 22 ) ) 
       return 0;    // ----------------------------------------------->


    // Group letters
    vGroupDef           aG; //(10);
    vColGroupDef        vColGroup1( vAli1[0].size(), -1 );
    vColGroupDef        vColGroup2( vAli2[0].size(), -1 );

    //+int    molType = f_getMolType( vAli1 ); // see above
    f_getColorGroups( cfg.sGroupFile, aG, molType );
    //getch();

    vGroupValueDef      vGroupValue1( aG.size(),0.);
    vGroupValueDef      vGroupValue2( aG.size(),0.);

    // Max group of columns and distribution of groups in entire align.
    f_statGroup( vAli1, aG, vColGroup1, vGroupValue1 ); 
    f_statGroup( vAli2, aG, vColGroup2, vGroupValue2 ); 

   oAli1.pGroup        = &aG;
   oAli2.pGroup        = &aG;


   oAli1.pAli2         = &vAli2;               // used only for comparison picture (simulate)
   oAli2.pAli2         = &vAli1;               // used only for comparison picture (simulate)

   oAli1.WeightLineSums2 =  sum2;
   oAli2.WeightLineSums2 =  sum1;

   oAli1.WeightLineMax2  =  oAli2.WeightLineMax;  
   oAli2.WeightLineMax2  =  oAli1.WeightLineMax;

   oAli1.WeightLineLen2  =  iSeqLen2;
   oAli2.WeightLineLen2  =  iSeqLen1;


   image_cfg.pColGroup     = &vColGroup1;
   image_cfg.pColGroup2    = &vColGroup2;
   image_cfg.pGroupValue   = &vGroupValue1;     // percent

   // cout << " -- weight max = " << oAli1.WeightLineMax << ", " << oAli2.WeightLineMax << endl;

   if( f_findType( cfg.vType, 21 ) ) {
       image_cfg.type = 21;

       image_cfg.pColGroup     = &vColGroup1;
       image_cfg.pColGroup2    = &vColGroup2;

       image_cfg.iH = iH1;
       image_cfg.WeightLineMax = oAli1.WeightLineMax;          // *** 22.08.22
       AligraImage( oAli1, image_cfg, (cfg.image_pref+string("ac1.png")).c_str(),
                (cfg.title).c_str(), (cfg.algo[0]).c_str() , (cfg.algo[1]).c_str() );

       image_cfg.pColGroup     = &vColGroup2;
       image_cfg.pColGroup2    = &vColGroup1;
       image_cfg.iH = iH2;
       image_cfg.WeightLineMax = oAli2.WeightLineMax;          // *** 22.08.22
       AligraImage( oAli2, image_cfg, (cfg.image_pref+string("bc1.png")).c_str(),
                (cfg.title).c_str(), (cfg.algo[1]).c_str() , (cfg.algo[0]).c_str() );
   }

   if( f_findType( cfg.vType, 22 ) ) {
       image_cfg.type = 22;

       image_cfg.pColGroup     = &vColGroup1;
       image_cfg.pColGroup2    = &vColGroup2;
       image_cfg.iH = iH1;
       image_cfg.WeightLineMax = oAli1.WeightLineMax;          // *** 22.08.22
       AligraImage( oAli1, image_cfg, (cfg.image_pref+string("ac2.png")).c_str(),
                (cfg.title).c_str(), (cfg.algo[0]).c_str() , (cfg.algo[1]).c_str() );

       image_cfg.pColGroup     = &vColGroup2;
       image_cfg.pColGroup2    = &vColGroup1;
       image_cfg.iH = iH2;
       image_cfg.WeightLineMax = oAli2.WeightLineMax;          // *** 22.08.22
       AligraImage( oAli2, image_cfg, (cfg.image_pref+string("bc2.png")).c_str(),
                (cfg.title).c_str(), (cfg.algo[1]).c_str() , (cfg.algo[0]).c_str() );
   }

   return 0;


/*****

   // Test values:  -------------------------------------------------


   vSNames.push_back( "Human" );              // 1
   vSNames.push_back( "Chicken" );            // 2
   vSNames.push_back( "Dogfish" );            // 3
   vSNames.push_back( "Lamprey" );            // 4
   vSNames.push_back( "Barley" );             // 5
   vSNames.push_back( "Maizey" );             // 6    'u'
   vSNames.push_back( "Lacto casei" );        // 7
   vSNames.push_back( "Bacillus stea" );      // 8
   vSNames.push_back( "Lacto plant" );        // 9
   vSNames.push_back( "Therma mari" );        // 10
   vSNames.push_back( "Bifido" );             // 11
   vSNames.push_back( "Thermus aqua" );       // 12
   vSNames.push_back( "Mycoplasma" );         // 13

   vSNames.push_back( "Lacto ca2"   );        // 14    'u'


   mOff[0]  = OffsetDescr( 'g' );
   mOff[1]  = OffsetDescr( '0' );
   mOff[80] = OffsetDescr( 'e' );   // 80 - just after the end

   vOffs.push_back( mOff );         // 1
   vOffs.push_back( mOff );         // 2
   vOffs.push_back( mOff );         // 3

   mOff.clear();
   mOff[-1] = OffsetDescr( 'm' );           // 'u'
   mOff[0]  = OffsetDescr( '0' );
   mOff[80] = OffsetDescr( 'e' );
   vOffs.push_back( mOff );         // 4

   mOff.clear();
   mOff[0]  = OffsetDescr( '0' );
   mOff[56] = OffsetDescr( 'd', -1 );
   mOff[57] = OffsetDescr( 'g' );
   mOff[58] = OffsetDescr( '0' );
   mOff[80] = OffsetDescr( 'e' );


   mOff[30] = OffsetDescr( 'd', -1 );
   mOff[31] = OffsetDescr( 'd', 2 );

   //mOff[32] = OffsetDescr( 'd', -10 );
   //mOff[33] = OffsetDescr( 'd', -50 );
   //mOff[34] = OffsetDescr( 'd', -5 );


   vOffs.push_back( mOff );         // 5

   mOff.clear();
   mOff[0]  = OffsetDescr( '0' );
   mOff[80] = OffsetDescr( 'e' );
   vOffs.push_back( mOff );         // 6

   mOff.clear();
   mOff[0]  = OffsetDescr( 'g' );
   mOff[1]  = OffsetDescr( '0' );
   mOff[50] = OffsetDescr( 'd', 1 );
   mOff[60] = OffsetDescr( 'g' );
   mOff[61] = OffsetDescr( '0' );
   mOff[79] = OffsetDescr( 'g' );
   mOff[80] = OffsetDescr( 'e' );
   vOffs.push_back( mOff );         // 7
   
   mOff.clear();
   mOff[0]  = OffsetDescr( 'g' );
   mOff[1]  = OffsetDescr( '0' );
   mOff[80] = OffsetDescr( 'e' );
   vOffs.push_back( mOff );          // 8
 
   mOff.clear();
   mOff[0]  = OffsetDescr( '0' );
   mOff[50] = OffsetDescr( 'd', -1 );
   mOff[60] = OffsetDescr( 'g' );
   mOff[61] = OffsetDescr( '0' );
   mOff[80] = OffsetDescr( 'e' );
   vOffs.push_back( mOff );         // 9
 
   mOff.clear();
   mOff[0]  = OffsetDescr( '0' );
   mOff[50] = OffsetDescr( 'd', -1 );
   mOff[60] = OffsetDescr( 'g' );
   mOff[61] = OffsetDescr( '0' );
   mOff[80] = OffsetDescr( 'e' );
   vOffs.push_back( mOff );         // 10
   
   mOff.clear();
   mOff[0]  = OffsetDescr( 'g' );
   mOff[1]  = OffsetDescr( '0' );
   mOff[80] = OffsetDescr( 'e' );
   vOffs.push_back( mOff );          // 11

   mOff.clear();
   mOff[0]  = OffsetDescr( '0' );
   mOff[61] = OffsetDescr( 'g' );
   mOff[62] = OffsetDescr( '0' );
   mOff[80] = OffsetDescr( 'e' );
   vOffs.push_back( mOff );          // 12
   
   mOff.clear();
   mOff[0]  = OffsetDescr( 'g' );
   mOff[1]  = OffsetDescr( '0' );
   mOff[80] = OffsetDescr( 'e' );
   vOffs.push_back( mOff );          // 13

   mOff.clear();
   mOff[-1] = OffsetDescr( 'u' );               // 'u'
   mOff[0]  = OffsetDescr( 'g' );
   mOff[1]  = OffsetDescr( '0' );
   mOff[60] = OffsetDescr( 'g' );
   mOff[61] = OffsetDescr( '0' );
   mOff[79] = OffsetDescr( 'g' );
   mOff[80] = OffsetDescr( 'e' );
   vOffs.push_back( mOff );         // 14 (7)



   //nph-malign:
   //sWL = ".*+.+**.*.**.+.+.+++.+++*+*+.*+*+.+++.+*+++*+.*+..*.....+.++.**....++++***+**+.*";
   //sWLleg=" .+*";

   //ClustalW:
   sWL = " :: ::* * ** :   :   :..* *  : *   .   ..  *: ..                    : .::: **  *";
   sWLleg=" .:*";

   oAli.pWeightLine = &sWL;
   oAli.pWLlegends  = &sWLleg;

   oAli.pSNames     = &vSNames;
   oAli.pOffs       = &vOffs;
 
   sTit = "" ; sAlgo1="2nd"; sAlgo2 = "1st";
   image_cfg.iH = 300; //250; // 150; 180; //300;   // 200  // 600
   //AlicompImage( vSNames, vOffs, image_cfg, "lact42.png", sTit, sAlgo1, sAlgo2 );
   AlicompImage( oAli, image_cfg, "lact48.png", sTit, sAlgo1, sAlgo2 );
   //AlicompImage( oAli, image_cfg, "lact46.png", "", "2nd", "1st " );


  return 0;
*****/
}