~ubuntu-branches/ubuntu/precise/exonerate/precise

/****************************************************************\

*                                                                *

*  C4 dynamic programming library - sub-alignment regions        *

*                                                                *

*  Guy St.C. Slater..   mailto:guy@ebi.ac.uk                     *

*  Copyright (C) 2000-2008.  All Rights Reserved.                *

*                                                                *

*  This source code is distributed under the terms of the        *

*  GNU General Public License, version 3. See the file COPYING   *

*  or http://www.gnu.org/licenses/gpl.txt for details            *

*                                                                *

*  If you use this code, please keep this notice intact.         *

*                                                                *

\****************************************************************/

#include "sar.h"

/**/

SAR_ArgumentSet *SAR_ArgumentSet_create(Argument *arg){

    register ArgumentSet *as;

    static SAR_ArgumentSet sas;

    if(arg){

        as = ArgumentSet_create("SAR Options");

        /**/

        ArgumentSet_add_option(as, '\0', "quality", "percent",

            "HSP quality threshold", "0",

            Argument_parse_float, &sas.hsp_quality);

        /**/

        Argument_absorb_ArgumentSet(arg, as);

        }

    return &sas;

    }

/**/

static gboolean SAR_region_valid_hsp_exit(Region *region, HSP *hsp){

#ifndef G_DISABLE_ASSERT

    register gint query_width = Region_query_end(region)

                              - hsp->query_start;

    register gint target_width = Region_target_end(region)

                               - hsp->target_start;

    register gint query_prefix = Region_query_end(region)

                               - hsp->query_start;

#endif /* G_DISABLE_ASSERT */

    /* Check region corner meets the HSP */

    g_assert(Region_is_valid(region));

    g_assert((query_width * HSP_target_advance(hsp))

          == (target_width * HSP_query_advance(hsp)));

    /* Check region corner is at HSP word boundary */

    g_assert(!(query_width % HSP_query_advance(hsp)));

    g_assert(!(target_width % HSP_target_advance(hsp)));

    /* Check region boundary is between cobs and HSP end */

    g_assert(query_prefix >= 0);

    g_assert(query_prefix <= HSP_query_cobs(hsp));

    return TRUE;

    }

static gboolean SAR_region_valid_hsp_entry(Region *region, HSP *hsp){

#ifndef G_DISABLE_ASSERT

    register gint query_width = HSP_query_end(hsp)

                              - region->query_start;

    register gint target_width = HSP_target_end(hsp)

                               - region->target_start;

    register gint query_suffix = HSP_query_end(hsp)

                               - region->query_start;

#endif /* G_DISABLE_ASSERT */

    g_assert(Region_is_valid(region));

    /* Check region corner meets the HSP */

    g_assert((query_width * HSP_target_advance(hsp))

          == (target_width * HSP_query_advance(hsp)));

    /* Check region corner is at HSP word boundary */

    g_assert(!(query_width % HSP_query_advance(hsp)));

    g_assert(!(target_width % HSP_target_advance(hsp)));

    /* Check region boundary is between HSP start and cobs */

    g_assert(query_suffix >= 0);

    g_assert(query_suffix <= (HSP_query_end(hsp)-HSP_query_cobs(hsp)));

    return TRUE;

    }

/**/

static gboolean SAR_Terminal_calculate_start_region(HSP *hsp,

            HPair *hpair, Heuristic_Range *range, C4_Scope scope,

            Region *region){

    register gint to_shrink;

    register gboolean at_query_edge = FALSE, at_target_edge = FALSE;

    Region outer_box;

    /* Find outer box (cobs -> corner) */

    outer_box.query_start = 0;

    outer_box.target_start = 0;

    outer_box.query_length = HSP_query_cobs(hsp);

    outer_box.target_length = HSP_target_cobs(hsp);

    /* Find inner_box (end point) */

    region->query_start = hsp->query_start;

    region->target_start = hsp->target_start;

    region->query_length = 0;

    region->target_length = 0;

    /* Grow inner_box by external range */

    region->query_start -= range->external_query;

    region->target_start -= range->external_target;

    region->query_length += range->external_query;

    region->target_length += range->external_target;

    /* Grow inner_box on HSPs by internal range */

    region->query_length += range->internal_query;

    region->target_length += range->internal_target;

    /* Trim inner_box to intersection (inner,outer) */

    if(region->query_start < outer_box.query_start){

        region->query_length -= (outer_box.query_start

                                 - region->query_start);

        region->query_start = outer_box.query_start;

        }

    if(region->target_start < outer_box.target_start){

        region->target_length -= (outer_box.target_start

                                  - region->target_start);

        region->target_start = outer_box.target_start;

        }

    to_shrink = Region_query_end(region)

              - Region_query_end(&outer_box);

    if(to_shrink > 0)

        region->query_length -= to_shrink;

    to_shrink = Region_target_end(region)

              - Region_target_end(&outer_box);

    if(to_shrink > 0)

        region->target_length -= to_shrink;

    /* Fail when region is empty */

    if(region->query_length <= 0)

        return FALSE;

    if(region->target_length <= 0)

        return FALSE;

    /* Check region is still on HSP boundary */

    g_assert(SAR_region_valid_hsp_exit(region, hsp));

    /* Check validity for non-local models */

    at_query_edge = (region->query_start == 0);

    at_target_edge = (region->target_start == 0);

    switch(scope){

        case C4_Scope_ANYWHERE:

            return TRUE;

            break;

        case C4_Scope_CORNER:

            if(at_query_edge && at_target_edge)

                return TRUE;

            break;

        case C4_Scope_EDGE:

            if(at_query_edge || at_target_edge)

                return TRUE;

            break;

        case C4_Scope_QUERY:

            if(at_query_edge)

                return TRUE;

            break;

        case C4_Scope_TARGET:

            if(at_target_edge)

                return TRUE;

            break;

        }

    return FALSE;

    }

static gboolean SAR_Terminal_calculate_end_region(HSP *hsp,

            HPair *hpair, Heuristic_Range *range, C4_Scope scope,

            Region *region){

    register gint to_shrink;

    register gboolean at_query_edge = FALSE, at_target_edge = FALSE;

    Region outer_box;

    /* Find outer box (cobs -> corner) */

    outer_box.query_start = HSP_query_cobs(hsp);

    outer_box.target_start = HSP_target_cobs(hsp);

    outer_box.query_length = hsp->hsp_set->query->len

                           - HSP_query_cobs(hsp);

    outer_box.target_length = hsp->hsp_set->target->len

                           - HSP_target_cobs(hsp);

    /**/

    /* Find inner_box (end point) */

    region->query_start = HSP_query_end(hsp);

    region->target_start = HSP_target_end(hsp);

    region->query_length = 0;

    region->target_length = 0;

    /* Grow inner_box by external range */

    region->query_length += range->external_query;

    region->target_length += range->external_target;

    /* Grow inner_box on HSPs by internal range */

    region->query_start -= range->internal_query;

    region->query_length += range->internal_query;

    region->target_start -= range->internal_target;

    region->target_length += range->internal_target;

    /* Trim inner_box to intersection (inner,outer) */

    if(Region_query_end(region) > Region_query_end(&outer_box)){

        region->query_length -= (Region_query_end(region)

                                - Region_query_end(&outer_box));

        }

    if(Region_target_end(region) > Region_target_end(&outer_box)){

        region->target_length -= (Region_target_end(region)

                                  - Region_target_end(&outer_box));

        }

    to_shrink = outer_box.query_start

              - region->query_start;

    if(to_shrink > 0){

        region->query_start += to_shrink;

        region->query_length -= to_shrink;

        }

    to_shrink = outer_box.target_start

              - region->target_start;

    if(to_shrink > 0){

        region->target_start += to_shrink;

        region->target_length -= to_shrink;

        }

    /* Fail when region is empty */

    if(region->query_length <= 0)

        return FALSE;

    if(region->target_length <= 0)

        return FALSE;

    /* Check region is still on HSP boundary */

    g_assert(SAR_region_valid_hsp_entry(region, hsp));

    /* Check validity for non-local models */

    at_query_edge = (Region_query_end(region)

                     == hsp->hsp_set->query->len);

    at_target_edge = (Region_target_end(region)

                     == hsp->hsp_set->target->len);

    switch(scope){

        case C4_Scope_ANYWHERE:

            return TRUE;

            break;

        case C4_Scope_CORNER:

            if(at_query_edge && at_target_edge)

                return TRUE;

            break;

        case C4_Scope_EDGE:

            if(at_query_edge || at_target_edge)

                return TRUE;

            break;

        case C4_Scope_QUERY:

            if(at_query_edge)

                return TRUE;

            break;

        case C4_Scope_TARGET:

            if(at_target_edge)

                return TRUE;

            break;

        }

    return FALSE;

    }

/**/

static C4_Score SAR_find_start_component(HPair *hpair,

                Region *region, HSP *hsp, C4_Calc *calc, gint *prefix){

    register C4_Score component = 0;

    register gint i, query_pos = hsp->query_start,

                     target_pos = hsp->target_start;

    register Region *calc_region;

    g_assert(prefix);

    (*prefix) = (Region_query_end(region) - hsp->query_start)

              / HSP_query_advance(hsp);

    g_assert((*prefix) >= 0);

    calc_region = Region_create(hsp->query_start, hsp->target_start,

                     (*prefix) * HSP_query_advance(hsp),

                     (*prefix) * HSP_target_advance(hsp));

    C4_Calc_init(calc, calc_region, hpair->user_data);

    for(i = 0; i < (*prefix); i++){

        component += C4_Calc_score(calc, query_pos, target_pos,

                                   hpair->user_data);

        query_pos += HSP_query_advance(hsp);

        target_pos += HSP_target_advance(hsp);

        }

    g_assert(component <= hsp->score);

    C4_Calc_exit(calc, calc_region, hpair->user_data);

    Region_destroy(calc_region);

    return component;

    }

static C4_Score SAR_find_end_component(HPair *hpair,

                Region *region, HSP *hsp, C4_Calc *calc, gint *suffix){

    register C4_Score component = 0;

    register gint i, query_pos, target_pos;

    register Region *calc_region;

    g_assert(suffix >= 0);

    g_assert(suffix);

    (*suffix) = (HSP_query_end(hsp) - region->query_start)

              / HSP_query_advance(hsp);

    calc_region = Region_create(region->query_start,

                                region->target_start,

                    (*suffix) * HSP_query_advance(hsp),

                    (*suffix) * HSP_target_advance(hsp));

    C4_Calc_init(calc, calc_region, hpair->user_data);

    query_pos = region->query_start;

    target_pos = region->target_start;

    for(i = 0; i < (*suffix); i++){

        component += C4_Calc_score(calc, query_pos, target_pos,

                                   hpair->user_data);

        query_pos += HSP_query_advance(hsp);

        target_pos += HSP_target_advance(hsp);

        }

    g_assert(component <= hsp->score);

    C4_Calc_exit(calc, calc_region, hpair->user_data);

    Region_destroy(calc_region);

    return component;

    }

/**/

static void SAR_HSP_quality(HPair *hpair, C4_Calc *calc,

                            HSP *hsp, C4_Score component,

                            gint start, gint length,

                            gint *half_score, gint *max_score){

    register gint i, query_pos, target_pos;

    g_assert(half_score);

    g_assert(max_score);

    query_pos = hsp->query_start + (start * HSP_query_advance(hsp));

    target_pos = hsp->target_start + (start * HSP_target_advance(hsp));

    (*half_score) = (*max_score) = 0;

    for(i = 0; i < length; i++){

        (*half_score) += HSP_get_score(hsp, query_pos, target_pos);

        (*max_score) += HSP_query_self(hsp, query_pos);

        query_pos += HSP_query_advance(hsp);

        target_pos += HSP_target_advance(hsp);

        }

    return;

    }

SAR_Terminal *SAR_Terminal_create(HSP *hsp, HPair *hpair,

                                  Heuristic_Match *match,

                                  gboolean is_start){

    register SAR_Terminal *sar_terminal;

    Region region;

    register C4_Score component;

    gint prefix, suffix;

    register SAR_ArgumentSet *sas = SAR_ArgumentSet_create(NULL);

    register gboolean is_valid = FALSE;

    gint half_score, max_score;

    register gint start, length;

    prefix = suffix = 0;

    if(is_start){

        is_valid = SAR_Terminal_calculate_start_region(hsp, hpair,

                       match->start_terminal->range,

                       hpair->heuristic->model->start_state->scope,

                       &region);

    } else { /* is end */

        is_valid = SAR_Terminal_calculate_end_region(hsp, hpair,

                       match->end_terminal->range,

                       hpair->heuristic->model->end_state->scope,

                       &region);

        }

    if(!is_valid)

        return NULL;

    if(is_start){

        component = SAR_find_start_component(hpair,

                                   &region, hsp, match->portal->calc,

                                   &prefix);

        start = prefix;

        length = hsp->cobs - prefix;

    } else {

        component = SAR_find_end_component(hpair,

                                   &region, hsp, match->portal->calc,

                                   &suffix);

        start = hsp->cobs;

        length = hsp->length - hsp->cobs - suffix;

        }

    if(length && (sas->hsp_quality > 0.0)){

        SAR_HSP_quality(hpair, match->portal->calc,

                        hsp, component, start, length,

                        &half_score, &max_score);

        if((((gdouble)half_score/(gdouble)max_score) * 100.0)

          < sas->hsp_quality)

            return NULL;

        }

    sar_terminal = g_new(SAR_Terminal, 1);

    sar_terminal->region = Region_copy(&region);

    sar_terminal->component = component;

    return sar_terminal;

    }

void SAR_Terminal_destroy(SAR_Terminal *sar_terminal){

    g_assert(sar_terminal);

    Region_destroy(sar_terminal->region);

    g_free(sar_terminal);

    return;

    }

C4_Score SAR_Terminal_find_bound(SAR_Terminal *sar_terminal,

                                 Heuristic_Bound *heuristic_bound){

    g_assert(sar_terminal->region->query_length >= 0);

    g_assert(sar_terminal->region->target_length >= 0);

    g_assert(sar_terminal->region->query_length

            <= heuristic_bound->query_range);

    g_assert(sar_terminal->region->target_length

            <= heuristic_bound->target_range);

    return heuristic_bound->matrix[sar_terminal->region->query_length]

                                  [sar_terminal->region->target_length]

           - sar_terminal->component;

    }

C4_Score SAR_Terminal_find_score(SAR_Terminal *sar_terminal,

                                 Optimal *optimal, HPair *hpair){

    return Optimal_find_score(optimal, sar_terminal->region,

                              hpair->user_data, hpair->subopt)

          - sar_terminal->component;

    }

/**/

static void SAR_reduce_mid_overlap(HPair *hpair, HSP *src, HSP *dst,

            Region *region, C4_Calc *src_calc, C4_Calc *dst_calc){

    register C4_Score src_total, dst_total, max_total;

    register gint src_query_pos, src_target_pos,

                  dst_query_pos, dst_target_pos;

    register gint max_src_query_pos, max_src_target_pos,

                  max_dst_query_pos, max_dst_target_pos,

                  max_dist_from_centre;

    /* Reverse up dst region, counting dst_total */

    if(!(region->query_length + region->target_length))

        return;

    src_total = dst_total = 0;

    dst_query_pos = Region_query_end(region)

                  - HSP_query_advance(dst);

    dst_target_pos = Region_target_end(region)

                   - HSP_target_advance(dst);

    while((dst_query_pos >= region->query_start)

       && (dst_target_pos >= region->target_start)

       && (dst_query_pos >= dst->query_start)

       && (dst_target_pos >= dst->target_start)){

        dst_total += C4_Calc_score(dst_calc,

                                   dst_query_pos, dst_target_pos,

                                   hpair->user_data);

        dst_query_pos -= HSP_query_advance(dst);

        dst_target_pos -= HSP_target_advance(dst);

        }

    dst_query_pos += HSP_query_advance(dst);

    dst_target_pos += HSP_target_advance(dst);

    /* Go down src region, storing max total */

    src_query_pos = region->query_start;

    src_target_pos = region->target_start;

    max_total = dst_total;

    max_src_query_pos = src_query_pos;

    max_src_target_pos = src_target_pos;

    max_dst_query_pos = dst_query_pos;

    max_dst_target_pos = dst_target_pos;

    max_dist_from_centre = Region_query_end(region)

                         - src_query_pos;

    while((src_query_pos < Region_query_end(region))

       && (src_target_pos < Region_target_end(region))

       && (src_query_pos < HSP_query_end(src))

       && (src_target_pos < HSP_target_end(src))){

        src_total += C4_Calc_score(src_calc,

                            src_query_pos, src_target_pos,

                            hpair->user_data);

        while((src_query_pos >= dst_query_pos)

           || (src_target_pos >= dst_target_pos)){

            dst_total -= C4_Calc_score(dst_calc,

                                       dst_query_pos, dst_target_pos,

                                       hpair->user_data);

            dst_query_pos += HSP_query_advance(dst);

            dst_target_pos += HSP_target_advance(dst);

            }

        if(max_total <= (src_total + dst_total)){

            /* If the score better

             * or the distance from the overlap centre is less

             */

            if((max_total < (src_total + dst_total))

            || (ABS(Region_query_end(region)-src_query_pos)

                < max_dist_from_centre)){

                max_dist_from_centre = ABS(Region_query_end(region)

                                           - src_query_pos);

                max_total = (src_total + dst_total);

                max_src_query_pos = src_query_pos;

                max_src_target_pos = src_target_pos;

                max_dst_query_pos = dst_query_pos;

                max_dst_target_pos = dst_target_pos;

                }

            }

        src_query_pos += HSP_query_advance(src);

        src_target_pos += HSP_target_advance(src);

        }

    /**/

    region->query_start = max_src_query_pos;

    region->target_start = max_src_target_pos;

    region->query_length = max_dst_query_pos

                         - max_src_query_pos;

    region->target_length = max_dst_target_pos

                          - max_src_target_pos;

    g_assert(Region_is_valid(region));

    return;

    }

static void SAR_find_end_box(HPair *hpair, HSP *src, HSP *dst,

                             Region *region, Region *cobs_box,

                             C4_Calc *src_calc, C4_Calc *dst_calc){

    register gint

        query_overlap = (HSP_query_end(src) - dst->query_start),

        target_overlap = (HSP_target_end(src) - dst->target_start);

    register gint src_query_move = 0, src_target_move = 0,

                  dst_query_move = 0, dst_target_move = 0;

    region->query_start = MIN(HSP_query_end(src), dst->query_start);

    region->target_start = MIN(HSP_target_end(src), dst->target_start);

    region->query_length = MAX(HSP_query_end(src), dst->query_start)

                         - region->query_start;

    region->target_length = MAX(HSP_target_end(src), dst->target_start)

                          - region->target_start;

    if((query_overlap > 0) || (target_overlap > 0)){

        /* Snap end_box to HSPs within cobs */

        src_query_move = region->query_start

                   - cobs_box->query_start;

        src_target_move = region->target_start

                    - cobs_box->target_start;

        if((src_query_move <= 0) || (src_target_move <= 0)){

            src_query_move = src_target_move = 0;

        } else {

            src_query_move -= (src_query_move % HSP_query_advance(src));

            src_target_move -= (src_target_move

                                % HSP_target_advance(src));

            if((src_query_move / HSP_query_advance(src))

             < (src_target_move / HSP_target_advance(src))){

                src_target_move = (src_query_move

                                   / HSP_query_advance(src))

                                * HSP_target_advance(src);

            } else {

                src_query_move = (src_target_move

                                  / HSP_target_advance(src))

                               * HSP_query_advance(src);

                }

            }

        /**/

        dst_query_move = Region_query_end(cobs_box)

                   - Region_query_end(region);

        dst_target_move = Region_target_end(cobs_box)

                    - Region_target_end(region);

        if((dst_query_move <= 0) || (dst_target_move <= 0)){

            dst_query_move = dst_target_move = 0;

        } else {

            dst_query_move -= (dst_query_move

                               % HSP_query_advance(dst));

            dst_target_move -= (dst_target_move

                                % HSP_target_advance(dst));

            if((dst_query_move / HSP_query_advance(dst))

             < (dst_target_move / HSP_target_advance(dst))){

                dst_target_move = (dst_query_move

                                   / HSP_query_advance(dst))

                            * HSP_target_advance(dst);

            } else {

                dst_query_move = (dst_target_move

                                  / HSP_target_advance(dst))

                               * HSP_query_advance(dst);

                }

            }

        region->query_start = cobs_box->query_start

                            + src_query_move;

        region->target_start = cobs_box->target_start

                             + src_target_move;

        region->query_length = Region_query_end(cobs_box)

                             - dst_query_move

                             - region->query_start;

        region->target_length = Region_target_end(cobs_box)

                             - dst_target_move

                             - region->target_start;

        g_assert(SAR_region_valid_hsp_entry(region, src));

        g_assert(SAR_region_valid_hsp_exit(region, dst));

        SAR_reduce_mid_overlap(hpair, src, dst, region,

                               src_calc, dst_calc);

        }

    g_assert(SAR_region_valid_hsp_entry(region, src));

    g_assert(SAR_region_valid_hsp_exit(region, dst));

    return;

    }

static gboolean SAR_find_cobs_box(HSP *src, HSP *dst, Region *region){

    region->query_start = HSP_query_cobs(src);

    region->target_start = HSP_target_cobs(src);

    region->query_length = HSP_query_cobs(dst)

                         - region->query_start;

    region->target_length = HSP_target_cobs(dst)

                          - region->target_start;

    if(region->query_length <= 0)

        return FALSE;

    if(region->target_length <= 0)

        return FALSE;

    return TRUE;

    }

static gboolean SAR_Join_calculate_region(HPair *hpair,

                                          HSP *src, HSP *dst,

                                          Region *region,

                                          Heuristic_Pair *pair){

    register gint to_shrink;

    Region outer_box;

    /* Find outer_box (cobs box) */

    if(!SAR_find_cobs_box(src, dst, &outer_box))

        return FALSE;

    /* Find inner_box (end_box) */

    SAR_find_end_box(hpair, src, dst, region, &outer_box,

                     pair->src->portal->calc, pair->dst->portal->calc);

    /* Check within external range */

    if(region->query_length > (pair->join->src_range->external_query

                             + pair->join->dst_range->external_query))

        return FALSE;

    if(region->target_length > (pair->join->src_range->external_target

                              + pair->join->dst_range->external_target))

        return FALSE;

    /* Grow inner_region by internal_range */

    region->query_start -= pair->join->src_range->internal_query;

    region->query_length += (pair->join->src_range->internal_query

                           + pair->join->dst_range->internal_query);

    region->target_start -= pair->join->src_range->internal_target;

    region->target_length += (pair->join->src_range->internal_target

                           + pair->join->dst_range->internal_target);

    /* Trim inner_box to intersection (inner,outer) */

    to_shrink = outer_box.query_start

              - region->query_start;

    if(to_shrink > 0){

        region->query_start += to_shrink;

        region->query_length -= to_shrink;

        }

    to_shrink = outer_box.target_start

              - region->target_start;

    if(to_shrink > 0){

        region->target_start += to_shrink;

        region->target_length -= to_shrink;

        }

    to_shrink = Region_query_end(region)

              - Region_query_end(&outer_box);

    if(to_shrink > 0)

        region->query_length -= to_shrink;

    to_shrink = Region_target_end(region)

              - Region_target_end(&outer_box);

    if(to_shrink > 0)

        region->target_length -= to_shrink;

    /* Check region not empty */

    if(region->query_length < 1)

        return FALSE;

    if(region->target_length < 1)

        return FALSE;

    /* Check region is still on HSP boundaries */

    g_assert(SAR_region_valid_hsp_entry(region, src));

    g_assert(SAR_region_valid_hsp_exit(region, dst));

    return TRUE;

    }

SAR_Join *SAR_Join_create(HSP *src_hsp, HSP *dst_hsp, HPair *hpair,

                          Heuristic_Pair *pair){

    register SAR_Join *sar_join;

    register C4_Score src_component, dst_component;

    register gint src_length, dst_length;

    gint prefix, suffix, src_half_score, dst_half_score,

                         src_max_score, dst_max_score;

    Region region;

    register SAR_ArgumentSet *sas = SAR_ArgumentSet_create(NULL);

    if(!SAR_Join_calculate_region(hpair, src_hsp, dst_hsp, &region,

                                  pair))

        return NULL;

    src_component = SAR_find_end_component(hpair,

              &region, src_hsp, pair->src->portal->calc, &suffix);

    dst_component = SAR_find_start_component(hpair,

              &region, dst_hsp, pair->dst->portal->calc, &prefix);

    src_length = src_hsp->length - src_hsp->cobs - suffix;

    dst_length = dst_hsp->cobs - prefix;

    if((src_length + dst_length) && (sas->hsp_quality > 0.0)){

        SAR_HSP_quality(hpair, pair->src->portal->calc,

                        src_hsp, src_component,

                        src_hsp->cobs, src_length,

                        &src_half_score, &src_max_score);

        SAR_HSP_quality(hpair, pair->dst->portal->calc,

                        dst_hsp, dst_component,

                        prefix, dst_length,

                        &dst_half_score, &dst_max_score);

        if(((((gdouble)(src_half_score+dst_half_score))

           / ((gdouble)(src_max_score+dst_max_score))) * 100.0)

          < sas->hsp_quality)

            return NULL;

        }

    sar_join = g_new(SAR_Join, 1);

    sar_join->region = Region_copy(&region);

    sar_join->src_component = src_component;

    sar_join->dst_component = dst_component;

    sar_join->pair = pair;

    return sar_join;

    }

void SAR_Join_destroy(SAR_Join *sar_join){

    g_assert(sar_join);

    Region_destroy(sar_join->region);

    g_free(sar_join);

    return;

    }

C4_Score SAR_Join_find_bound(SAR_Join *sar_join){

    g_assert(sar_join->region->query_length >= 0);

    g_assert(sar_join->region->target_length >= 0);

    g_assert(sar_join->region->query_length

            <= sar_join->pair->join->bound->query_range);

    g_assert(sar_join->region->target_length

            <= sar_join->pair->join->bound->target_range);

    return sar_join->pair->join->bound->matrix