src/common/rotamer.cc

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <assert.h>
#include <string>

#include "open_prospect.h"
#include "tree_decomp.h"
#include "sm_matrix.h"



int ResidueAtoms::GetAtomIndex( int atom_name_in ) {
        int cur_atom = -1;
        for ( int i = 0; i < atom_count; i++ ) {
                if ( atom_name[i] == atom_name_in ) 
                        cur_atom = i;
        }       
        return cur_atom;        
}

void ResidueAtoms::SetAtomLoc( int atom_name_in, vec3 loc) {
        
        int cur_atom = -1;
        for ( int i = 0; i < atom_count; i++ ) {
                if ( atom_name[i] == atom_name_in ) 
                        cur_atom = i;
        }       
        if ( cur_atom == -1 ) { 
                cur_atom = atom_count;
                atom_count++;
                if ( atom_name_in >= 0 && atom_name_in < 4 ) {
                        bb_num[ atom_name_in ] = cur_atom;
                }
                atom_name[ cur_atom ] = atom_name_in;
        }
        
        atom[ cur_atom ][ 0 ] = loc[0];
        atom[ cur_atom ][ 1 ] = loc[1];
        atom[ cur_atom ][ 2 ] = loc[2];
}

void ResidueAtoms::SetTempFactor( int atom_index, float new_temp) {
        if (atom_index < atom_count) {
                temp_factor[ atom_index ] = new_temp;
        }
}


void ResidueAtoms::ReadAtomLine( char * line ) {
        char tmp_x_str[30], tmp_y_str[30], tmp_z_str[30];
        strcpy(tmp_x_str, "");
        strcpy(tmp_y_str, "");
        strcpy(tmp_z_str, "");
        
        if ( !strncmp(line, "ATOM", 4) || !strncmp(line, "HETATM", 6) ) {                       
                char tmp_str[50];
                strncpy( tmp_str, &(line[22]), 5 );
                tmp_str[5] = 0;
                
                if ( atom_count < kMaxResAtoms) {
                        //load in the atom num
                        strncpy( tmp_str, &(line[6]), 4 );
                        atom_num[ atom_count ] = atoi( tmp_str );
                        
                        strncpy(tmp_x_str, &(line[30]), 8);
                        tmp_x_str[8] = 0;
                        strncpy(tmp_y_str, &(line[38]), 8);
                        tmp_y_str[8] = 0;
                        strncpy(tmp_z_str, &(line[46]), 8);
                        tmp_z_str[8] = 0;
                        //location information/////////////////////////////
                        atom[atom_count][0] = atof(tmp_x_str);
                        atom[atom_count][1] = atof(tmp_y_str);
                        atom[atom_count][2] = atof(tmp_z_str);
                        
                        //Temperature factor/////////////////////////////
                        strncpy( tmp_str, &(line[60]), 6);
                        tmp_str[6] = 0;
                        temp_factor[ atom_count ] = atof( tmp_str );
                        
                        //residue id//////////////////////////////////////
                        strncpy(tmp_str, &(line[22]), 4 );
                        tmp_str[4] = 0;
                        res_num = atoi( tmp_str );
                        res_mod = line[26];
                        
                        //the_data->atom[i].residue[5] = 0;
                        //clean_space( the_data->atom[i].residue );
                        
                        //atom name//////////////////////////////////////
                        int atom_namenum = PDBAtomtoNum( &(line[13]) );
                        if ( atom_namenum >= 0 ) {
                                atom_name[ atom_count ] = atom_namenum;
                                if ( atom_namenum == anN || atom_namenum == anCa || 
                                         atom_namenum == anC || atom_namenum == anO  ) {
                                        bb_num[ atom_namenum ] = atom_count;
                                }
                                atom_count++;
                        }
                }
        }       
}




char *ResidueAtoms::WriteAtomLine(int i, char RS, int res_num, char chain) {
        if ( i > atom_count )
                return NULL;
        char *buffer = (char *)malloc( 200 );
        char tmp[10];
        AA1toAA3(RS, tmp);                      
        if (chain == 0)
                chain = ' ';
        sprintf(buffer, "ATOM  %5d  %-3s %+3s %c%5d   %8.3f%8.3f%8.3f  1.00%6.2f\n", 
                        atom_num[i], 
                        pdb_atomname[ atom_name[i] ], 
                        tmp, 
                        chain, 
                        res_num,
                        atom[i][0], 
                        atom[i][1], 
                        atom[i][2],
                        temp_factor[i] );
        return buffer;
}




void smTranslateCoordinates(vec3 pre_anchor[3], vec3 cur_anchor[3], vec3 A[], int N) {
        
        int i;
        vec3 pre_position, cur_position;
        float matrix1[3][3], matrix2[3][3], matrix3[3][3], matrix4[3][3];
        
        smCopy(pre_anchor[anCa], pre_position);
        smCopy(cur_anchor[anCa], cur_position);
        smMinus(cur_position);
        
        //      smComputeConvertMatrix(pre_anchor[0], pre_anchor[1], pre_anchor[2], matrix1);
        //      smComputeConvertMatrix(cur_anchor[0], cur_anchor[1], cur_anchor[2], matrix2);
        
        smComputeConvertMatrix(pre_anchor[anCa], pre_anchor[anN], pre_anchor[anC], matrix1);
        smComputeConvertMatrix(cur_anchor[anCa], cur_anchor[anN], cur_anchor[anC], matrix2);
        
        
        smInverseMatrix(matrix2, matrix3);
        smMultMat(matrix3, matrix1, matrix4);
        //      smMultiMatrix(matrix2, matrix1, matrix4, 3);
        
        for (i=0; i<N; i++)
                smAdd(cur_position, A[i]);
        for (i=0; i<N; i++)
                smMultVec(A[i], matrix4);
        for (i=0; i<N; i++)
                smAdd(pre_position, A[i]);
        return;
}


int Prospect_Target_Rotamer(ProspectParam *param, ResidueAtoms *structure, char residue, int rotamer_number) {
        //FIXME: need to check this code for compatibility with new residue structure format    
        //first step, copy in residue  lib co-ordinates
        vec3 a[20];
        int res_num = AA1toAANum(residue);
        if (res_num >= 0 && res_num < 20) {
                ResidueInfoElement *res_info = &param->residue_info[ res_num ];                 
                //copy atoms not part of the backbone
                int atom_count = 0;
                for (int i = 0; i < res_info->structure.atom_count; i++) {
                        if ( res_info->structure.atom_name[i] > anBBcount ) {
                                a[ atom_count ][0] = res_info->structure.atom[i][0];
                                a[ atom_count ][1] = res_info->structure.atom[i][1];
                                a[ atom_count ][2] = res_info->structure.atom[i][2];                    
                                atom_count++;
                        }
                }               
                assert( atom_count < 20);
                
                //apply rotamer angle adjustments               
                
                //translate the rotamer to the backbone
                smTranslateCoordinates( structure->atom, res_info->structure.atom, a, atom_count);              
                for (int i = 0; i < atom_count; i++) {
                        structure->atom[i ][0] = a[i][0]; //FIXME: need offset here
                        structure->atom[i ][1] = a[i][1];
                        structure->atom[i ][2] = a[i][2];
                }               
                if ( residue != 'G' && residue != 'P' ) {
                        RotamerElement *rotamer = &param->rotamer_array[ res_num ].rotamer[ rotamer_number ];
                        //      assert( res_info->torsion_count == param->rotamer_array[ res_num ].chi_count );
                        for (int i = 0; i < res_info->torsion_count; i++) {
                                //                      for (int i = 0; i < res_info->torsion_count && i < 1; i++) 
                                int tor_a = res_info->torsion_a[i];
                                int tor_b = res_info->torsion_b[i];
                                
                                vec3 axis;
                                smSub(structure->atom[tor_b], structure->atom[tor_a], axis);
                                smNormalize(axis);
                                float rot_mat[3][3];
                                smMakeRotationMatrix( (rotamer->chi[i] - res_info->torsion_val[i])/180.0*M_PI, axis, rot_mat);
                                assert( !isnan(rot_mat[0][0]) );
                                //      smMakeRotationMatrix(  res_info->torsion_val[i] - rotamer->chi[i], axis, rot_mat);
                                /* Do rotation */
                                for(int j=0; j < structure->atom_count; j++) {
                                        if ( res_info->torsion_mask[i][j] ) {
                                                vec3 tmp;
                                                tmp[0] = structure->atom[j][0] - structure->atom[tor_a][0];
                                                tmp[1] = structure->atom[j][1] - structure->atom[tor_a][1];
                                                tmp[2] = structure->atom[j][2] - structure->atom[tor_a][2];
                                                structure->atom[j][0] = rot_mat[0][0]*tmp[0] + rot_mat[0][1]*tmp[1] + rot_mat[0][2]*tmp[2] + structure->atom[tor_a][0];
                                                structure->atom[j][1] = rot_mat[1][0]*tmp[0] + rot_mat[1][1]*tmp[1] + rot_mat[1][2]*tmp[2] + structure->atom[tor_a][1];
                                                structure->atom[j][2] = rot_mat[2][0]*tmp[0] + rot_mat[2][1]*tmp[1] + rot_mat[2][2]*tmp[2] + structure->atom[tor_a][2];
                                                assert( !isnan(structure->atom[j][0]) && !isnan(structure->atom[j][1]) && !isnan(structure->atom[j][2]));
                                        }
                                }
                        }
                }               
        }
        return 0;
}


pValType calc_single_energy(TargetStruct *target_data, ResidueAtoms *res_struct) {
        pValType energy = 0.0;
        for (int atom_a = 0;
                 atom_a < res_struct->atom_count; atom_a++) {
                if ( res_struct->atom_name[ atom_a ] >= anBBcount ) {
                        for (int res_b = 0; res_b < target_data->len; res_b++) {
                                for (int atom_b = 0; atom_b < target_data->structure[res_b].atom_count; atom_b++) {
                                        if ( res_struct->atom_name[ atom_b ] < anBBcount ) {
                                                energy += calc_LRET( res_struct->atom[atom_a], 
                                                                                         target_data->structure[ res_b ].atom[atom_b],
                                                                                         res_struct->radius[atom_a], 
                                                                                         target_data->structure[ res_b ].radius[atom_b] );
                                        }
                                }
                        }
                }
        }
        return energy;
}



pValType calc_pair_energy(ResidueAtoms *res_a, ResidueAtoms *res_b) {
        pValType energy = 0.0;
        for (int atom_a = 0;
                 atom_a < res_a->atom_count; atom_a++) {
                if (res_a->atom_name[ atom_a ] >= anBBcount) {
                        for (int atom_b = 0;
                                 atom_b < res_b->atom_count; atom_b++) {
                                if (res_a->atom_name[ atom_b ] >= anBBcount) {
                                        energy += calc_LRET( res_a->atom[atom_a], 
                                                                                 res_b->atom[atom_b], 
                                                                                 res_a->radius[atom_a], 
                                                                                 res_b->radius[atom_b] );
                                }
                        }
                }
        }
        return energy;
}



//#define target_hash( x, y )  (( x * target_data->len + y))




int rot_mask_active(RotamerMask *mask, int i) {
        int count = 0;
        for (int j = 0; j < mask->count[i]; j++) {
                if ( mask->state[i][j] ) 
                        count++;
        }
        return count;
}




int Prospect_Target_RotamerMask_Init(ProspectParam *param, TargetStruct *target_data, RotamerMask *rotamer_mask) {
        int total_count = 0;
        
        rotamer_mask->count = (long *)malloc(sizeof(long) * target_data->len);
        for (int i = 0; i < target_data->len; i++) {
                rotamer_mask->count[i] = 0;
        }
        if (  target_data->structure == NULL ) 
                return 1;
        
        for (int i = 0; i < target_data->len; i++) {
                if ( target_data->structure[i].atom_count >= 4 ) {
                        char res_num = AA1toAANum(target_data->residue[i].RS) ;
                        for (int j = 0; j < param->rotamer_array[ res_num ].rotamer_count; j++) {
                                if ( fabs(target_data->residue[i].psi - param->rotamer_array[ res_num ].rotamer[j].psi) < param->rotamer_array[res_num].psi_bin &&
                                         fabs(target_data->residue[i].phi - param->rotamer_array[ res_num ].rotamer[j].phi) < param->rotamer_array[res_num].phi_bin ) {
                                        rotamer_mask->count[i]++;
                                        total_count++;
                                }
                        }
                } else {
                        rotamer_mask->count[i] = 0;             
                }
        }
        rotamer_mask->len = target_data->len;
        rotamer_mask->state = (char **)malloc(sizeof(char *) * target_data->len);
        rotamer_mask->state[0] = (char *)malloc(sizeof(char) * total_count);
        rotamer_mask->rot_num = (long **)malloc(sizeof(long *) * target_data->len);
        rotamer_mask->rot_num[0] = (long *)malloc(sizeof(long) * total_count);
        rotamer_mask->structure = (ResidueAtoms **)malloc(sizeof(ResidueAtoms *) * target_data->len);
        rotamer_mask->structure[0] = (ResidueAtoms *)malloc(sizeof(ResidueAtoms) * total_count);
        rotamer_mask->single_score = (pValType **)malloc(sizeof(pValType *) * target_data->len);
        rotamer_mask->single_score[0] = (pValType *)malloc(sizeof(pValType) * total_count);
        
        int cur_pos = 0;
        for (int i = 0; i < target_data->len; i++) {
                if ( target_data->structure[i].atom_count >= 4) {
                        rotamer_mask->state[i] = rotamer_mask->state[0] + cur_pos;
                        rotamer_mask->structure[i] = rotamer_mask->structure[0] + cur_pos;
                        rotamer_mask->single_score[i] = rotamer_mask->single_score[0] + cur_pos;
                        rotamer_mask->rot_num[i] = rotamer_mask->rot_num[0] + cur_pos;
                }
                cur_pos += rotamer_mask->count[i];
        }
        
        for (int i = 0; i < target_data->len; i++) {
                if  ( target_data->structure[i].atom_count >= 4) {
                        
                        char res_num = AA1toAANum(target_data->residue[i].RS) ;
                        int rot_count = 0;
                        for (int j = 0; j < param->rotamer_array[ res_num ].rotamer_count; j++) {
                                if ( fabs(target_data->residue[i].psi - param->rotamer_array[ res_num ].rotamer[j].psi) < param->rotamer_array[res_num].psi_bin &&
                                         fabs(target_data->residue[i].phi - param->rotamer_array[ res_num ].rotamer[j].phi) < param->rotamer_array[res_num].phi_bin ) {
                                        rotamer_mask->rot_num[i][ rot_count ] = j;
                                        rot_count++;
                                }
                        }
                        for (int j =0; j < rotamer_mask->count[i]; j++) {
                                rotamer_mask->state[i][j] = 1;
                                rotamer_mask->structure[i][j] = target_data->structure[i];
                                Prospect_Target_Rotamer(param, &rotamer_mask->structure[i][j], target_data->residue[i].RS, rotamer_mask->rot_num[i][j]);
                                rotamer_mask->single_score[i][j] = calc_single_energy(target_data, &rotamer_mask->structure[i][j]);
                        }
                }
        }
        
        
        
        rotamer_mask->pair_score = (pValType ***)malloc(sizeof(pValType **) * target_data->len * target_data->len);
        for (int i = 0; i < target_data->len; i++) {
                if  ( target_data->structure[i].atom_count >= 4)  {
                        for (int j = 0; j < target_data->len; j++) {
                                if  ( j > i ) {
                                        if (  ( target_data->structure[i].atom_count >= 4)  && i != j) {
                                                //int offset =  target_hash( i, j);
                                                int offset = i * target_data->len + j;
                                                rotamer_mask->pair_score[offset] = (pValType **)malloc(sizeof(pValType *) * rotamer_mask->count[i]  );
                                                rotamer_mask->pair_score[offset][0] = (pValType *)malloc(sizeof(pValType) * rotamer_mask->count[i] * rotamer_mask->count[j]);
                                                for (int  k = 0; k < rotamer_mask->count[i]; k++ ) {
                                                        rotamer_mask->pair_score[offset][k] = rotamer_mask->pair_score[offset][0] + k * rotamer_mask->count[j];  
                                                }
                                        } else {
                                                //rotamer_mask->pair_score[ target_hash( i, j) ] = NULL;                                        
                                                int offset = i * target_data->len + j;
                                                rotamer_mask->pair_score[ offset ] = NULL;                      
                                        }
                                } else {
                                        int offset = i * target_data->len + j;
                                        rotamer_mask->pair_score[ offset ] = NULL;
                                }
                        }
                } else {
                        for (int j = 0; j < target_data->len; j++) {
                                //rotamer_mask->pair_score[ target_hash( i, j) ] = NULL;
                                int offset = i * target_data->len + j;
                                rotamer_mask->pair_score[ offset ] = NULL;
                        }
                }
        }               
        
        for (int i = 0; i < target_data->len; i++) {
                if  ( target_data->structure[i].atom_count >= 4)  {
                        for (int j = i+1; j < target_data->len; j++) {
                                if  ( target_data->structure[i].atom_count >= 4)  {
                                        float dist2 = pow( target_data->structure[i].atom[ anCb ][0] - target_data->structure[j].atom[ anCb ][0], 2) +
                                        pow( target_data->structure[i].atom[ anCb ][1] - target_data->structure[j].atom[ anCb ][1], 2) +
                                        pow( target_data->structure[i].atom[ anCb ][2] - target_data->structure[j].atom[ anCb ][2], 2) ;
                                        if ( dist2 < 240.25) {  //15.5^2  angstoms
                                                for (int i_r = 0; i_r < rotamer_mask->count[i]; i_r++) {
                                                        //int offset = target_hash( i, j);
                                                        for (int j_r = 0; j_r < rotamer_mask->count[j]; j_r++) {
                                                                //rotamer_mask->pair_score[offset][i_r][j_r] = 
                                                                //calc_pair_energy(&rotamer_mask->structure[i][i_r], &rotamer_mask->structure[j][j_r]);
                                                                rot_mask_pairscore_set( rotamer_mask, i, j, i_r, j_r, 
                                                                                                                calc_pair_energy(&rotamer_mask->structure[i][i_r],
                                                                                                                                                 &rotamer_mask->structure[j][j_r]) 
                                                                                                                );
                                                        }
                                                }               
                                        } else {
                                                for (int i_r = 0; i_r < rotamer_mask->count[i]; i_r++) {
                                                        for (int j_r = 0; j_r < rotamer_mask->count[j]; j_r++) {
                                                                rot_mask_pairscore_set( rotamer_mask, i, j, i_r, j_r, 0.0 );
                                                        }
                                                }
                                        }
                                }               
                        }                       
                }
        }       
        
        
        
        /*
         for (int i = 0; i < target_data->len; i++) {
                 for (int j =0; j < rotamer_mask->count[i]; j++) {
                         rotamer_mask->pair_score[i][j][
                                 */     
        

        return 0;
}


void Prospect_Target_RotamerMask_Free(RotamerMask *mask) {
        free( mask->count );
        free(mask->state[0] );
        free(mask->state );
        free(mask->rot_num[0] );
        free(mask->rot_num  );
        free(mask->structure[0] );
        free(mask->structure  );
        free(mask->single_score[0] );
        free(mask->single_score  );
        for (int i = 0; i < mask->len; i++) {
                for (int j = 0; j < mask->len; j++) {
                        if ( mask->pair_score[ i * mask->len + j ] != NULL )
                                free( mask->pair_score[ i * mask->len + j ] );
                }
        }
}

pValType rot_mask_pairscore_get( RotamerMask *mask, long vertex_1, long vertex_2, long state_1, long state_2) {
        if (vertex_2 > vertex_1) {
                return mask->pair_score[ vertex_1 * mask->len + vertex_2 ][ state_1 ][ state_2 ];               
        }
        return mask->pair_score[ vertex_2 * mask->len + vertex_1 ][ state_2 ][ state_1 ];       
}


void rot_mask_pairscore_set( RotamerMask *mask, long vertex_1, long vertex_2, long state_1, long state_2, pValType score) {
        if (vertex_2 > vertex_1) {
                mask->pair_score[ vertex_1 * mask->len + vertex_2 ][ state_1 ][ state_2 ] = score;              
        } else {
                mask->pair_score[ vertex_2 * mask->len + vertex_1 ][ state_2 ][ state_1 ] = score;      
        }
}


/*
 void Prospect_Target_RotamerMask_Compress(RotamerMask *mask) {
         for (int i = 0; i < mask->len; i++) {
                 int k = 0; 
                 for (int j = 0; j < mask->count[i]; j++) {
                         if ( mask->state[i][j] ) {
                                 mask->state[i][k] = mask->state[i][j];
                                 mask->rot_num[i][k] = mask->rot_num[i][j];
                                 mask->structure[i][k] = mask->structure[i][j];
                                 k++;
                         }
                 }
         }
 }
 */


int Prospect_RotamerMask_DEE(RotamerMask *mask) {
        
                
        char removal_found = 0; 
        do {
                removal_found = 0;
                for (int i = 0; i < mask->len; i++) {
                        if ( mask->count[i] != 0) {
                                //fprintf(stderr, "Residue %d (%d rotamers)\n", i, mask->count[i]);
                                for (int rot_a = 0; rot_a < mask->count[i]; rot_a++) {
                                        if ( mask->state[i][rot_a] ) {
                                                float e_single_a = mask->single_score[i][rot_a];
                                                for (int rot_b = 0; rot_b < mask->count[i]; rot_b++) {
                                                        if (rot_a != rot_b && mask->state[i][rot_b]) {
                                                                float e_single_b = mask->single_score[i][rot_b];                                                        
                                                                float e_pair_sum = 0;
                                                                for ( int res_c = 0; res_c < mask->len; res_c++ ) {
                                                                        if ( res_c != i && mask->count[res_c] != 0 ) {
                                                                                float min = 1000000000;
                                                                                //int offset = i * mask->len + res_c;//target_hash( i, res_c );
                                                                                for ( int rot_c = 0; rot_c < mask->count[res_c]; rot_c++) {
                                                                                        if ( mask->state[res_c][rot_c] ) {                                                                                                                                      
                                                                                                //float e_pair_a = mask->pair_score[ offset ][ rot_a ][ rot_c ];
                                                                                                //float e_pair_b = mask->pair_score[ offset ][ rot_b ][ rot_c ];
                                                                                                float e_pair_a = rot_mask_pairscore_get( mask, i, res_c, rot_a, rot_c );
                                                                                                float e_pair_b = rot_mask_pairscore_get( mask, i, res_c, rot_b, rot_c );                                                                                                                                                                                                                                                                         
                                                                                                float tmp = e_pair_b - e_pair_a;
                                                                                                if ( tmp < min ) {
                                                                                                        min = tmp;
                                                                                                }
                                                                                        }
                                                                                }
                                                                                e_pair_sum += min;
                                                                        }
                                                                }                                                       
                                                                if ( e_single_b - e_single_a + e_pair_sum > 0 ) {
                                                                        //fprintf(stderr, "Eliminate res: %d rotamer: %d by %d (%f)\n", i, rot_b, rot_a, e_single_b - e_single_a + e_pair_sum);
                                                                        mask->state[i][rot_b] = 0;
                                                                        removal_found = 1;
                                                                }
                                                        }
                                                }
                                        }
                                }
                        } else {
                                //fprintf(stderr, "Residue %d skipped\n", i);
                        }
                }
        } while (removal_found);        
        return 0;
}


/*
int Prospect_RotamerMask_DEE2(RotamerMask *mask) {
        
        char removal_found = 0;
        
        do {
                removal_found = 0;
                for (int i = 0; i < mask->len; i++) {
                        int i_count = rot_mask_active(mask, i);                 
                        if ( i_count > 1 ) {
                                for (int j = 0; j < mask->len; j++) {
                                        int j_count = rot_mask_active(mask, j);
                                        if (i != j && j_count > 1) {
                                                for ( int ir_1 = 0; ir_1 < mask->count[i]; ir_1++) {
                                                        if ( mask->state[i][ ir_1 ] ) {
                                                                for ( int jr_1 = 0; jr_1 < mask->count[j]; jr_1++) {
                                                                        if ( ir_1 != jr_1 && mask->state[j][ jr_1 ] ) {
                                                                                long offset_1 = i * mask->len + j;
                                                                                float val_1 = mask->single_score[i][ir_1] + mask->single_score[j][jr_1] + mask->pair_score[offset_1][ir_1][jr_1];
                                                                                for ( int ir_2 = 0; ir_2 < mask->count[i]; ir_2++) {
                                                                                        if ( ir_1 != ir_2 && mask->state[i][ ir_2 ] ) {
                                                                                                for ( int jr_2 = 0; jr_2 < mask->count[j]; jr_2++) {
                                                                                                        if ( jr_1 != jr_2 && mask->state[j][ jr_2 ] ) {
                                                                                                                float val_2 = mask->single_score[i][ir_2] + mask->single_score[j][jr_2] + mask->pair_score[offset_1][ir_2][jr_2];
                                                                                                                
                                                                                                                for (int k = 0; k < mask->len; k++) {
                                                                                                                        if (k != i && k != j ) {
                                                                                                                                
                                                                                                                        }
                                                                                                                }       
                                                                                                        }
                                                                                                }
                                                                                        }
                                                                                }
                                                                        }
                                                                }
                                                        }
                                                }
                                        }
                                }
                        }
                }
        } while (removal_found);
        return 0;
}
*/

int Prospect_RotamerMask_DEE_Split(RotamerMask *mask) {
        
        int max_rotcount = 0;
        for (int i = 0; i < mask->len; i++) {
                if ( max_rotcount < mask->count[i])
                        max_rotcount = mask->count[i];
        }
        
        float y_val[ mask->len][ max_rotcount ];                
        char y_val_set[ mask->len][ max_rotcount ];     

        char removal_found = 0; 
        do {
                removal_found = 0;
                for (int i = 0; i < mask->len; i++) {
                        for ( int ir_1 = 0; ir_1 < mask->count[i]; ir_1++) {
                                if ( mask->state[i][ ir_1 ] ) {
                                        for ( int ir_2 = 0; ir_2 < mask->count[i]; ir_2++) {
                                                if ( ir_1 != ir_2 && mask->state[i][ ir_2 ] ) {
                                                        for (int j = 0; j < mask->len; j++) {
                                                                if ( i != j ) {
                                                                        y_val[j][ir_2] = 1e20;
                                                                        y_val_set[j][ir_2] = 0;
                                                                        for ( int jr_1 = 0; jr_1 < mask->count[j]; jr_1++) {
                                                                                if ( mask->state[j][ jr_1 ] ) {
                                                                                        float val1 = rot_mask_pairscore_get( mask, i, j, ir_1, jr_1 );
                                                                                        float val2 = rot_mask_pairscore_get( mask, i, j, ir_2, jr_1 );
                                                                                        float test_val = val1 - val2;
                                                                                        if ( test_val < y_val[j][ir_2] ) {
                                                                                                y_val[j][ir_2] = test_val;
                                                                                                y_val_set[j][ir_2] = 1;
                                                                                        }
                                                                                }
                                                                        }
                                                                }
                                                        }                                                       
                                                }
                                        }                               
                                        for (int k = 0; k < mask->len; k++) {
                                                if ( k != i && rot_mask_active(mask, k) > 0 ) {
                                                        char v_mask[ max_rotcount ];
                                                        for (int tmp = 0; tmp < mask->count[k]; tmp++) 
                                                                v_mask[tmp] = 0;
                                                        for ( int ir_2 = 0; ir_2 < mask->count[i]; ir_2++) {
                                                                if ( ir_1 != ir_2 && mask->state[i][ ir_2 ] ) {
                                                                        float val_x = mask->single_score[i][ir_1] - mask->single_score[i][ir_2];
                                                                        for (int j = 0; j < mask->len; j++) {
                                                                                if ( j != i && j != k) {
                                                                                        if ( y_val_set[j][ir_2] )
                                                                                                val_x += y_val[j][ir_2];
                                                                                }
                                                                        }
                                                                        for ( int kr_1 = 0; kr_1 < mask->count[k]; kr_1++) {
                                                                                if ( mask->state[k][ kr_1 ] ) {
                                                                                        float val1 = rot_mask_pairscore_get( mask, i, k, ir_1, kr_1 );
                                                                                        float val2 = rot_mask_pairscore_get( mask, i, k, ir_2, kr_1 );
                                                                                        if ( val_x + (val1 - val2) > 0 ) {
                                                                                                v_mask[ kr_1 ] = 1;
                                                                                        }
                                                                                }
                                                                        }
                                                                }
                                                        }
                                                        char eliminate = 1;
                                                        for (int tmp = 0; tmp < mask->count[k]; tmp++) {
                                                                if ( mask->state[k][tmp] && !v_mask[ tmp ]) {
                                                                        eliminate = 0;
                                                                }
                                                        }
                                                        if (eliminate) {
                                                                mask->state[i][ir_1] = 0;
                                                                k = mask->len;
                                                                removal_found = 1;
                                                        }
                                                }
                                        }
                                }                       
                        }
                }
        } while (removal_found);        
        return 0;
}


void sc_tree_print(char **matrix, long matrix_size, long cur_node, long source) {
        int i;
        
        for (i = 0; i < matrix_size; i++) {
                if (cur_node != i && i != source && matrix[cur_node][i]) {
                        fprintf(stderr, "%ld -> %d\n", cur_node, i);
                        sc_tree_print(matrix, matrix_size, i, cur_node);                        
                }
        }
}


int Prospect_RotamerMask_TreeDecomp( RotamerMask *mask ) {      
        
        //figure out the tree decomposition of the problem      
        //first, create the graph of connected points
        long vertex_count = 0;
        for (int i = 0; i < mask->len; i++) {
                if ( mask->count[i] != 0) {
                        int count = 0;
                        for (int rot_a = 0; rot_a < mask->count[i]; rot_a++) {
                                if ( mask->state[i][rot_a] ) {
                                        count++;
                                }
                        }
                        char connect = 0;
                        for (int j = 0; j < mask->len; j++) {
                                if ( j != i ) {
                                        //int offset = i * mask->len + j;//target_hash( i, j );
                                        for (int ir =0; ir < mask->count[i]; ir++) {
                                                if ( mask->state[i][ir])  {
                                                        for (int jr =0; jr < mask->count[j]; jr++) {
                                                                if ( mask->state[j][jr])  {
                                                                        //if (mask->pair_score[ offset ][ ir ][ jr ] > 0 ) 
                                                                        if ( rot_mask_pairscore_get( mask, i, j, ir, jr ) > 0 ) {
                                                                                connect = 1;
                                                                        }
                                                                }
                                                        }
                                                }
                                        }
                                }
                        }
                        if (connect) {
                                vertex_count++;                         
                        }
                }
        }
        
        
        //do the vertex hashes for the residues, and count the number of states
        long *vertex_hash = (long *)malloc(sizeof(long) * vertex_count );
        long *state_count = (long *)malloc(sizeof(long) * vertex_count );
        long total_state_count = 0;
        vertex_count = 0;       
        for (int i = 0; i < mask->len; i++) {
                if ( mask->count[i] != 0) {
                        int count = 0;
                        for (int rot_a = 0; rot_a < mask->count[i]; rot_a++) {
                                if ( mask->state[i][rot_a] ) {
                                        count++;
                                }
                        }
                        char connect = 0;
                        for (int j = 0; j < mask->len; j++) {
                                if ( j != i ) {
//                                      int offset = i * mask->len + j;//target_hash( i, j );
                                        for (int ir =0; ir < mask->count[i]; ir++) {
                                                if ( mask->state[i][ir])  {
                                                        for (int jr =0; jr < mask->count[j]; jr++) {
                                                                if ( mask->state[j][jr])  {
                                                                        //if (mask->pair_score[ offset ][ ir ][ jr ] > 0 ) 
                                                                        if (rot_mask_pairscore_get( mask, i, j, ir, jr) > 0 ) {
                                                                                connect = 1;
                                                                        }
                                                                }
                                                        }
                                                }
                                        }
                                }
                        }
                        if (connect) {
                                vertex_hash[ vertex_count ] = i;
                                state_count[ vertex_count ] = count;
                                total_state_count += count;
                                vertex_count++;
                        }
                }
        }
        
        //we also need to create the rotamer hash
        long **rot_hash = (long **)malloc(sizeof(long *) * vertex_count);
        rot_hash[0] = (long *)malloc(sizeof(long) * total_state_count);
        total_state_count = 0;
        for (int i = 0; i < vertex_count; i++) {
                rot_hash[i] = rot_hash[0] + total_state_count;
                int k = 0;
                for (int j = 0; j <  mask->count[ vertex_hash[i] ]; j++) {
                        if ( mask->state[ vertex_hash[i] ][j] ) {
                                rot_hash[i][k] = j;
                                k++;
                        }
                }
                total_state_count += state_count[ i ];
        }
        
        
        //draw the graph connection matrix
        char **graph_matrix = td_matrix_init( vertex_count );
        for (int i = 0; i < vertex_count; i++) {
                for (int j = i+1; j < vertex_count; j++) {
                        int vi = vertex_hash[i];
                        int vj = vertex_hash[j];                        
                        //int offset = vi * mask->len + vj;   //target_hash( vi, vj )
                        char connect = 0;
                        for (int ir = 0; ir < state_count[i]; ir++) {
                                int vir = rot_hash[i][ir];
                                for (int jr =0; jr < state_count[j]; jr++) {
                                        int vjr = rot_hash[j][jr];
                                        //if (mask->pair_score[ offset ][ vir ][ vjr ] > 0 ) 
                                        if ( rot_mask_pairscore_get( mask, vi, vj, vir, vjr ) > 0 ) 
                                                connect = 1;
                                }
                        }
                        if (connect) {
                                graph_matrix[i][j] = 1;
                                graph_matrix[j][i] = 1;
                        }
                }
        }       
        
        td_graph *sc_graph = td_graph_init( vertex_count, state_count, graph_matrix );  
        
        //fill in the score values of the graph
        for (int i = 0; i < vertex_count; i++) {
                int vi = vertex_hash[i];
                for (int ir = 0; ir < state_count[i]; ir++) {
                        int vir = rot_hash[i][ir];
                        sc_graph->vertex_score[i][ir] = mask->single_score[ vi ][vir];
                }               
                for (int j = i+1; j < vertex_count; j++) {
                        if ( sc_graph->matrix[i][j]) {
                                int vj = vertex_hash[j];
                                long graph_offset = td_edge_offset( sc_graph, i, j );
                                //long mask_offset = vi * mask->len + vj;  //target_hash( vi, vj );
                                for (int ir = 0; ir < state_count[i]; ir++) {
                                        int vir = rot_hash[i][ir];
                                        for (int jr =0; jr < state_count[j]; jr++) {
                                                int vjr = rot_hash[j][jr];
                                                //sc_graph->edge_score[ graph_offset ][ ir ][ jr ] = mask->pair_score[ mask_offset ][ vir ][ vjr ];
                                                sc_graph->edge_score[ graph_offset ][ ir ][ jr ] = rot_mask_pairscore_get( mask, vi, vj, vir, vjr ) ;
                                        }
                                }
                        }
                }
        }
        
        
        
        td_decomp* decomp = td_graph_decomp_tri( sc_graph );
        
#ifdef DEBUG
        for (int i = 0; i < decomp->bag_count; i++) {
                fprintf(stderr, "bag %d: ", i);
                for (int j = 0; j < decomp->bag_size[i]; j++) {
                        fprintf(stderr, "%d ", decomp->bag[i][j] );
                }
                fprintf(stderr, "\n");
        }
        fprintf(stderr, "Tree root: %d\n", decomp->root);
        sc_tree_print(decomp->bag_matrix, decomp->bag_count, decomp->root, -1);
#endif
        
#if 1
        /*
         for (int i = 0; i < decomp->bag_count; i++) {
                 printf("%d: ", i);
                 for ( int j = 0; j < decomp->bag_size[i]; j++) {
                         printf(" %d", vertex_hash[ decomp->bag[i][j] ]);
                 }
                 printf("\n");
         }
         */
        long *opt_state = td_graph_decomp_solve( decomp );
        
        for (int i = 0; i < vertex_count; i++) {
                int vi = vertex_hash[i];
                int vir = rot_hash[i][ opt_state[i] ];
                for (int j = 0; j < mask->count[vi]; j++) {
                        if ( j == vir ) {
                                mask->state[vi][j] = 1;
                        } else {
                                mask->state[vi][j] = 0;
                        }                               
                }
                
        }
        
#else
        
        //output the graph
        printf("graph G {\n");
        /*
         for (int i = 0; i < mask->len; i++) {
                 if ( mask->count[i] != 0) {
                         int count = 0;
                         int a = 0;
                         for (int rot_a = 0; rot_a < mask->count[i]; rot_a++) {
                                 if ( mask->state[i][rot_a] ) {
                                         if (count == 0)
                                                 a = rot_a;
                                         count++;
                                 }
                         }
                         //     if ( count > 1 ) {
                         for (int j = i+1; j < mask->len; j++) {
                                 if ( j != i ) {
                                         int offset = i * mask->len + j;//target_hash( i, j );
                                         char connect = 0;
                                         for (int ir =0; ir < mask->count[i]; ir++) {
                                                 if ( mask->state[i][ir])  {
                                                         for (int jr =0; jr < mask->count[j]; jr++) {
                                                                 if ( mask->state[j][jr])  {
                                                                         if (mask->pair_score[ offset ][ ir ][ jr ] > 0 ) {
                                                                                 connect = 1;
                                                                         }
                                                                 }
                                                         }
                                                 }
                                         }
                                         if (connect) {
                                                 printf("\t%d -- %d;\n", i, j);
                                         }
                                 }
                                 //     }
                         }
                 }
         }              
         for (int i = 0; i < mask->len; i++) {
                 if ( mask->count[i] != 0) {
                         int count = 0;
                         int a = 0;
                         for (int rot_a = 0; rot_a < mask->count[i]; rot_a++) {
                                 if ( mask->state[i][rot_a] ) {
                                         if (count == 0)
                                                 a = rot_a;
                                         count++;
                                 }
                         }
                         if (count == 1 || count == 0) {
                                 printf("\t%d [shape=triangle];\n", i);
                         }
                 }
         }
         */
        
        for (int i = 0; i < vertex_count; i++) {
                for (int j = i+1; j < vertex_count; j++) {
                        if (graph_matrix[i][j]) {
                                printf("\t%d -- %d\n", i, j);
                        }
                }               
        }
        
        /*      
                for (int i = 0; i < decomp->bag_count; i++) {
                        printf("\tsubgraph subgraph%d {\n", i);
                        for (int j = 0; j < decomp->bag_size[i]; j++) {
                                for ( int k = j+1; k < decomp->bag_size[i]; k++) {
                                        if ( graph_matrix[ decomp->bag[i][j] ][ decomp->bag[i][k] ] ) {
                                                printf( "\t\tsub%d_%d -- sub%d_%d;\n", i, decomp->bag[i][j], i, decomp->bag[i][k] );
                                        }
                                }
                        }
                        printf("\t}\n");
                }       
         */
        for (int i = 0; i < decomp->bag_count; i++) {
                for (int j = 0; j < decomp->bag_count; j++) {
                        if ( decomp->bag_matrix[i][j] ) {
                                printf("\tbag%d -- bag%d;\n", i, j);
                        }
                }
        }
        for (int i = 0; i < decomp->bag_count; i++) {
                char buffer[2000];
                strcpy(buffer, "");
                for (int j = 0; j < decomp->bag_size[i]; j++) {
                        char buffer2[100];
                        sprintf(buffer2, "%d, ", decomp->bag[i][j]);
                        strcat(buffer, buffer2);
                }
                printf("\tbag%d [label=\"%s\"];\n", i, buffer);
        }
        
        printf("}\n");  
        
        
#endif
        return 0;
}

---