interval_set.cc

Go to the documentation of this file.

// Interval_set Class -*- C++ -*-

// $Id: interval_set.cc 927 2009-11-12 17:10:57Z mihaly $
// Copyright (C) 2005 Ferenc Domes
//
// This file is part of the COCONUT API.  This library
// is free software; you can redistribute it and/or modify it under the
// terms of the Library GNU General Public License as published by the
// Free Software Foundation; either version 2, or (at your option)
// any later version.

// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// Library GNU General Public License for more details.

// As a special exception, you may use this file as part of a free software
// library without restriction.  Specifically, if other files instantiate
// templates or use macros or inline functions from this file, or you compile
// this file and link it with other files to produce an executable, this
// file does not by itself cause the resulting executable to be covered by
// the Library GNU General Public License.  This exception does not however
// invalidate any other reasons why the executable file might be covered by
// the Library GNU General Public License.

//------------------------------------------------------------------------------
// IMPLEMENTAION
//------------------------------------------------------------------------------

#include <coconut_config.h>
#include <coconut_vect.h>
#include <interval_set.h>
// #include <map>
#include <stdlib.h>
#include <cstring>

namespace coco {

static int sortint(const void *a, const void *b)
{
  if(((const interval*)a)->sup() < ((const interval*)b)->inf())
    return -1;
  else if(((const interval*)a)->inf() > ((const interval*)b)->sup())
    return 1;
  else
    return 0;
}

void interval_set::sort(){
        if (length>1)
                qsort((void*)data, length, sizeof(interval), sortint);
}

interval_set& interval_set::operator=(const interval_set &a){
        if(this != &a){
                length=a.length;
                maxlength=a.maxlength;extension=a.extension;
                sorting=a.sorting;
                if(data) delete [] data;
                data = new interval[maxlength];
                for(int i=0;i<length;i++)
                        data[i]=a.data[i];
        }
        return (*this);
}

void interval_set::setSorting (bool s){
        if (!sorting && s && length>1)
                sort();
        sorting=s;
}

bool interval_set::is_subset(const interval_set &a) const{
        for(int i=0;i<a.getLength();i++){
                bool subs=false;
                for(int j=0;j<length;j++){
                        if (data[j].subset(a[i])){
                                 subs=true;break;
                        }
                }
                if (!subs) return false;
                
        }
        return true;
}

bool interval_set::is_subset(const interval &a) const{
        return is_subset(interval_set(a));
}

void interval_set::add(const interval &a)
{
        int ainfin(-1), asupin(-1), i, j;
        for(i=0;i<length;i++){
                if (data[i].contains(a))
                        return;
                else if (a.contains(data[i]))
                        remove(i--);
                else if (data[i].contains(a.inf()))
                        ainfin=i;
                else if (data[i].contains(a.sup())){
                    if(sorting)
                    {
                        if(ainfin > 0)
                            data[ainfin]=interval(data[ainfin].inf(),
                                data[i].sup());
                        else
                            data[i]=interval(a.inf(),data[i].sup());
                        return;
                    }
                    else
                        asupin=i;
                }
                else if(sorting && ainfin > 0 && data[i].inf() > a.sup())
                {
                    data[ainfin]=interval(data[ainfin].inf(), a.sup());
                    return;
                }
        }
        if(ainfin > 0)
        {
          if(asupin > 0)
          {
            data[ainfin]=interval(data[ainfin].inf(),data[asupin].sup());
            remove(asupin);
          }
          else
            data[ainfin]=interval(data[ainfin].inf(),a.sup());
          return;
        }
        else if(asupin > 0)
        {
          data[asupin]=interval(a.inf(), data[asupin].sup());
          return;
        }
        if (length==maxlength){
                interval *data2=new interval[maxlength+extension];
                for (i=0;i<maxlength;i++)
                        data2[i]=data[i];
                if(data) delete[] data;
                data=data2;
                maxlength+=extension;
        }
        int pos=length;
        if (sorting){
                pos=0;
                while(pos<length){
                        if(a < data[pos])
                                break;
                        pos++;
                }
                if (pos<length && length>0){
                        for(int pcp=length-1; pcp>pos; --pcp)
                          data[pcp]=data[pcp-1];
                }
        }
        data[pos]=a;
        if((unsigned int)length == MAXLEN)
        {
                if(!sorting) sort(); // otherwise, this does not work
                typedef std::pair<double,int> pair_double_int;
                __vectgen(pair_double_int,gaps,fillextension);
                for(i=0; i<(int)fillextension; ++i)
                  gaps[i].first = -COCO_INF;
                for(i=0; i<length; ++i)
                {
                  double thisgap = data[i+1].inf()-data[i].sup();
                  for(j=0; j<(int)fillextension; ++j)
                    if(thisgap > gaps[j].first)
                      break;
                  if(j<(int)fillextension)
                  {
                    for(int k=fillextension-1; k>j; --k)
                      gaps[k]=gaps[k-1];
                    gaps[j]=std::make_pair(thisgap,i);
                  }
                }
                __vectgen(int,removeit,length);
                memset(removeit, 0, length*sizeof(int));
                for(i=0; i<(int)fillextension; ++i)
                {
                  if(gaps[i].first == -COCO_INF)
                    break;
                  removeit[gaps[i].second]=1;
                }
                length -= i;
                for(i=0, j=0; i<length; ++i, ++j)
                {
                  if(removeit[j])
                  {
                    data[i]=interval(data[j].inf(),data[j+1].sup());
                    ++j;
                  }
                  else if(i != j)
                    data[i]=data[j];
                }
        }
        length++;
}

interval interval_set::remove (int pos)
{
        if(pos<0 || pos>=length)
                return interval(NAN,NAN);
        interval r(data[pos]);
        if (length==maxlength-extension && length > 0)
        {       // realloc to smaller size
                maxlength-=extension;
                interval *data2=new interval[maxlength];
                for (int i=0;i<pos;i++){
                        data2[i]=data[i];
                }
                for (int i=pos+1;i<length;i++){
                        data2[i-1]=data[i];
                }
                if(data) delete[] data;
                data=data2;
        }
        else    // just remove the entry;
        {
                for(int i=pos+1; i<length; ++i)
                        data[i-1]=data[i];
        }
        --length;
        return r;
}

interval_set interval_set::remove (const interval &a){
        interval_set rem;
        for (int i=0;i<length;i++){
                if(a.contains(data[i])){
                        rem.add(data[i]);
                        remove(i--);
                }
                else if (data[i].contains(a)){
                        rem.add(a);
                        if(data[i].inf() < a.inf())
                        {
                                data[i]=interval(data[i].inf(),a.inf());
                                if(data[i].sup() > a.sup())
                                        add(interval(a.sup(),data[i].sup()));
                        }
                        else // data[i].sup() > a.sup()
                                data[i]=interval(a.sup(),data[i].sup());
                        return rem;
                }
                else if (data[i].inf() < a.inf() && data[i].sup() > a.inf()){
                        rem.add(interval(a.inf(),data[i].sup()));
                        data[i]=interval(data[i].inf(),a.inf());
                }
                else if (data[i].inf() < a.sup() && data[i].sup() > a.sup()){
                        rem.add(interval(data[i].inf(),a.sup()));
                        data[i]=interval(a.sup(),data[i].sup());
                }
        }
        return rem;
}

interval interval_set::minimum() const{
        if(length==0)
                return 0;
        else if (sorting)
                return data[0];
        else{
                interval m(COCO_INF);
                for(int i=0;i<length;i++){
                        if (data[i]<m)
                                m=data[i];
                }
                return m;
        }
}

interval interval_set::maximum() const{
        if(length==0)
                return 0;
        else if (sorting)
                return data[length-1];
        else{
                interval m(-COCO_INF);
                for(int i=0;i<length;i++){
                        if (m<data[i])
                                m=data[i];
                }
                return m;
        }
}

interval interval_set::hull() const{
        if(length==0)
                return interval::EMPTY();
        else
                return interval(minimum().inf(),maximum().sup());
}

double interval_set::width() const{
        return hull().width();
}

interval_set& interval_set::round_to_integer()
{
  int i, j;
  for (i=0, j=-1; i<length; ++i)
  {
    data[i].round_to_integer();
    if(!data[i].is_empty())
    {
      ++j;
      if(i != j)
        data[j] = data[i];
    }
  }
  length = j+1;
  return *this;
}

interval_set& interval_set::intersectwith(const interval_set &a){
        interval_set i2;
        for (int i=0;i<length;i++){
                for (int j=0;j<a.getLength();j++){
                        interval I(data[i]);
                        I.intersectwith(a[j]);
                        if (!I.empty())
                                i2.add(I);
                }
        }
        return *this=i2;
}

void interval_set::intersect_power(const interval_set &a, int n){
        interval_set i2;
        for (int i=0;i<length;i++){
                for (int j=0;j<a.getLength();j++){
                        interval I(data[i]);
                        I.intersect_power(a[j],n);
                        if (!I.empty())
                                i2.add(I);
                }
        }
        *this=i2;
}

void interval_set::intersect_div(const interval_set &a, const interval_set &b){
        interval_set i2;
        for (int j=0;j<a.getLength();j++){
                for (int k=0;k<b.getLength();k++){
                        interval_set J(divide(a[j],b[k]));
                        J.intersectwith(*this);
                        if (!J.is_empty())
                                i2.add(J);
                }
        }
        *this=i2;
}

void interval_set::intersect_invsqr_wc(const interval_set &a, double b, double c){
        interval_set i2;
        for (int i=0;i<length;i++){
                for (int j=0;j<a.getLength();j++){
                        interval_set up, lo;
                        up = interval_set((sqrt(a[j]) - c)/b);
                        lo = interval_set((-sqrt(a[j]) - c)/b);
                        up.intersectwith(*this);
                        lo.intersectwith(*this);
                        if (!up.is_empty())
                                i2.add(up);
                        if (!lo.is_empty())
                                i2.add(lo);
                }
        }
        *this=i2;
}

void interval_set::intersect_invpower_wc(const interval_set &a, double b, int n){
        interval_set i2;
        for (int i=0;i<length;i++){
                for (int j=0;j<a.getLength();j++){
                        interval I(data[i]);
                        I.intersect_invpower_wc(a[j],b,n);
                        if (!I.empty())
                                i2.add(I);
                }
        }
        *this=i2;
}

void interval_set::intersect_invsin_wc(const interval_set &a, double b, double c){
        interval_set i2;
        for (int i=0;i<length;i++){
                for (int j=0;j<a.getLength();j++){
                        interval I(data[i]);
                        I.intersect_invsin_wc(a[j],b,c);
                        if (!I.empty())
                                i2.add(I);
                }
        }
        *this=i2;
}

void interval_set::intersect_invcos_wc(const interval_set &a, double b, double c){
        interval_set i2;
        for (int i=0;i<length;i++){
                for (int j=0;j<a.getLength();j++){
                        interval I(data[i]);
                        I.intersect_invcos_wc(a[j],b,c);
                        if (!I.empty())
                                i2.add(I);
                }
        }
        *this=i2;
}

void interval_set::intersect_invgauss_wc(const interval_set &a, double b, double c, double d){
        interval_set i2;
        for (int i=0;i<length;i++){
                for (int j=0;j<a.getLength();j++){
                        interval I(data[i]);
                        I.intersect_invgauss_wc(a[j],b,c,d);
                        if (!I.empty())
                                i2.add(I);
                }
        }
        *this=i2;
}

interval_set& interval_set::operator +=(const interval_set &a){
        *this=*this+a;
        return *this;
}

interval_set& interval_set::operator -=(const interval_set &a){
        *this=*this-a;
        return *this;
}

interval_set& interval_set::operator *=(const interval_set &a){
        *this=*this*a;
        return *this;
}

interval_set& interval_set::operator /=(const interval_set &a){
        *this=*this/a;
        return *this;
}

double width(const interval_set& a){return a.width();}
//OPERATORS

interval_set divide(const interval &a,const interval &b)
{
        interval_set c;
        bool x2_used;
        c.add(division_part1(a, b, x2_used));
        if(x2_used){
                c.add(division_part2(a, b));
        }
        return c;
}

bool operator== (const interval_set &a,const interval_set &b){
        for(int i=0;i<a.getLength();i++){
                for(int j=0;j<b.getLength();j++){
                        if (a[i] != b[j])
                                return false;
                }
        }
        return true;    
}

//FUNCTIONS
interval_set imin(const interval_set &a, const interval_set &b){
        interval_set i2;
        for (int i=0;i<a.getLength();i++){
                for (int j=0;j<b.getLength();j++){
                        interval I=imin(a[i],b[j]);
                        if (!I.empty())
                                i2.add(I);
                }
        }
        return i2;
}
interval_set imax(const interval_set &a, const interval_set &b){
        interval_set i2;
        for (int i=0;i<a.getLength();i++){
                for (int j=0;j<b.getLength();j++){
                        interval I=imax(a[i],b[j]);
                        if (!I.empty())
                                i2.add(I);
                }
        }
        return i2;
}
interval_set power(const interval_set &a, int n){
        interval_set i2;
        for (int i=0;i<a.getLength();i++){
                interval I=power(a[i],n);
                if (!I.empty())
                        i2.add(I);
        }
        return i2;
}

interval_set sqr(const interval_set &a){
        interval_set i2;
        for (int i=0;i<a.getLength();i++){
                interval I=sqr(a[i]);
                if (!I.empty())
                        i2.add(I);
        }
        return i2;
}
interval_set sqrt(const interval_set &a){
        interval_set i2;
        for (int i=0;i<a.getLength();i++){
                interval I=sqrt(a[i]);
                if (!I.empty())
                        i2.add(I);
        }
        return i2;
}
interval_set exp(const interval_set &a){
        interval_set i2;
        for (int i=0;i<a.getLength();i++){
                interval I=exp(a[i]);
                if (!I.empty())
                        i2.add(I);
        }
        return i2;
}
interval_set log(const interval_set &a){
        interval_set i2;
        for (int i=0;i<a.getLength();i++){
                interval I=log(a[i]);
                if (!I.empty())
                        i2.add(I);
        }
        return i2;
}
interval_set abs(const interval_set &a){
        interval_set i2;
        for (int i=0;i<a.getLength();i++){
                interval I=abs(a[i]);
                if (!I.empty())
                        i2.add(I);
        }
        return i2;
}
interval_set sin(const interval_set &a){
        interval_set i2;
        for (int i=0;i<a.getLength();i++){
                interval I=sin(a[i]);
                if (!I.empty())
                        i2.add(I);
        }
        return i2;
}
interval_set cos(const interval_set &a){
        interval_set i2;
        for (int i=0;i<a.getLength();i++){
                interval I=cos(a[i]);
                if (!I.empty())
                        i2.add(I);
        }
        return i2;
}
interval_set pow(const interval_set &a, const interval_set &b){
        interval_set i2;
        for (int i=0;i<a.getLength();i++){
                for (int j=0;j<b.getLength();j++){
                        interval I=pow(a[i],b[j]);
                        if (!I.empty())
                                i2.add(I);
                }
        }
        return i2;
}
interval_set atan2(const interval_set &a, const interval_set &b){
        interval_set i2;
        for (int i=0;i<a.getLength();i++){
                for (int j=0;j<b.getLength();j++){
                        interval I=atan2(a[i],b[j]);
                        if (!I.empty())
                                i2.add(I);
                }
        }
        return i2;
}

interval_set round_to_integer(const interval_set &a){
        interval_set i2(a);
        i2.round_to_integer();
        return i2;
}

} // end of namespace coco

namespace std {
ostream& operator<< (ostream &s, const coco::interval_set &a){
        if (a.getLength() == 0)
                s<< "{ empty }";
        else{
                bool first(true);
                s<< "{";
                for (int i=0;i<a.getLength();i++)
                {
                  if(first) first=false; else s << ",";
                        s << a[i];
                }
                s<< "}";
        }
        return s;
}

} // end of namespace std

---