massTool.cpp

Go to the documentation of this file.

/* 
 * Copyright (C) 2017 D Haley
 * 
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 * 
 * This program 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
 * GNU General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */

#include "atomprobe/algorithm/massTool.h"

#include <cmath>
#include <limits>
#include <algorithm>
#include <numeric>
#include <utility>
#include <map>


#include "atomprobe/helper/misc.h"
#include "atomprobe/helper/aptAssert.h"

using std::vector;
using std::make_pair;
using std::map;
using std::pair;

using namespace AtomProbe;

//Holder that allows us to track position in weight stack when solving
// knapsack problem
class WEIGHT_SEARCH_ENTRY
{
    public:
    vector<Weight> curWeights,allowableWeights;
    size_t offset;
    float cumulativeWeight;
    
};

//Simple template class to create a stack of known maximum size 
template<class T>
class FixedStack
{
    T *array;
    size_t curSize;
public: 
    FixedStack(size_t size);
    FixedStack(const FixedStack &f); //Not implemented
    ~FixedStack();
    
    T &top();
    void pop();
    void push(const T &t);

    bool empty() const {return !curSize;}
};

template<class T>
FixedStack<T>::FixedStack(size_t size)
{
    array = new T[size];
    curSize=0;
}

template<class T>
void FixedStack<T>::push(const T &t)
{
    curSize++;
    array[curSize]=t;
}

template<class T>
void FixedStack<T>::pop()
{
    ASSERT(curSize);
    curSize--;
}

template<class T>
T& FixedStack<T>::top()
{
    return array[curSize];
}

template<class T>
FixedStack<T>::~FixedStack()
{
    delete[] array;
}

void MassTool::preprocessKnapsack(vector<Weight> &weights, float totalWeight, float tolerance, unsigned int maxCombine)
{

    //Delete weights that cannot possibly be part of the solution
    // * Weight is over the target + tolerance
    vector<bool> badWeight(weights.size(),false);
    for(size_t ui=0;ui<weights.size();ui++)
    {
        if(weights[ui].mass > totalWeight+tolerance)
            badWeight[ui]=true;
    }
    
    
    vectorMultiErase(weights,badWeight);

}

void MassTool::bruteKnapsack(vector<Weight> &weights, float targetWeight, 
        float tolerance, unsigned int maxObjects, vector<vector<Weight> >  &solutions)
{
    using std::pair;


    if(!weights.size())
        return;

    //Reject weights that are not part of the solution, to reduce 
    // the search space
    preprocessKnapsack(weights,targetWeight,tolerance,maxObjects);
    
    //Loop over all weights, adding objects until we are within |tolerant weight + targetweight|
    const float upperWeight = targetWeight+tolerance;
    const float lowerWeight = targetWeight-tolerance;

    std::sort(weights.begin(),weights.end());
    
    WEIGHT_SEARCH_ENTRY weightEntry;
    weightEntry.offset=0;
    weightEntry.cumulativeWeight=0;
    weightEntry.allowableWeights=weights;
    
    FixedStack<WEIGHT_SEARCH_ENTRY> searchStack(maxObjects+1);
    searchStack.push(weightEntry);
    
    //---

    while(!searchStack.empty())
    {   
        vector<Weight> allowableWeights,curWeights;
        
        unsigned insertOffset;

        //Check the top of the stack
        //--
        insertOffset=searchStack.top().offset;

    
        //Don't proceed if we have used all the child weights   
        if(insertOffset >=searchStack.top().allowableWeights.size())
        {
            //We have run out of allowable weights to use
            searchStack.pop();
            continue;
        }

        //Tell current stack to move to next child
        searchStack.top().offset++;

        //Copy the top of the stack
        curWeights=searchStack.top().curWeights;
        allowableWeights=searchStack.top().allowableWeights;
        
        //Insert a weight into the sack
        float remainingWeight,currentWeight;
        curWeights.push_back(allowableWeights[insertOffset]);
        currentWeight=searchStack.top().cumulativeWeight +allowableWeights[insertOffset].mass;

        if( currentWeight <=upperWeight &&  currentWeight >=lowerWeight)
        {
            //TODO: Use solution tree, to more efficiently store possible solutions?
            solutions.push_back(curWeights);
        }
            

        //Erase disallowed weights.
        //the weights are sorted in increasing size, so just
        // cut once we hit the right threshold
        //--------
        remainingWeight=upperWeight-currentWeight;
        if(remainingWeight > 0)
        {
            for(size_t ui=0;ui<allowableWeights.size();ui++)
            {
                if(allowableWeights[ui].mass > remainingWeight)
                {
                    allowableWeights.resize(ui+1);
                    break;
                }
            }
        }   
        else 
            allowableWeights.clear();
        //--------

        //Add the next weight to the stack, if we have not run out of number of objects
        if(allowableWeights.size() && curWeights.size() != maxObjects)
        {
            weightEntry.curWeights.swap(curWeights);
            weightEntry.allowableWeights.swap(allowableWeights);
            weightEntry.cumulativeWeight=currentWeight;
            //Search in decreasing weight only, to not generate duplicates
            weightEntry.offset=insertOffset;
            searchStack.push(weightEntry);
        }


    }

}

void MassTool::bruteKnapsack(vector<Weight> &weights, const vector<float> &targetWeight, 
        float tolerance, unsigned int maxObjects, vector<vector<Weight> >  &solutions)
{
    using std::pair;


    if(weights.empty() || targetWeight.empty())
        return;

    float maxWeight=*(std::max_element(targetWeight.begin(),targetWeight.end()));

    //Reject weights that are not part of the solution, to reduce 
    // the search space
    preprocessKnapsack(weights,maxWeight,tolerance,maxObjects);
    
    //Loop over all weights, adding objects until we are within |tolerant weight + targetweight|
    const float upperWeight = maxWeight+tolerance;

    std::sort(weights.begin(),weights.end());
    
    WEIGHT_SEARCH_ENTRY weightEntry;
    weightEntry.offset=0;
    weightEntry.cumulativeWeight=0;
    weightEntry.allowableWeights=weights;
    
    FixedStack<WEIGHT_SEARCH_ENTRY> searchStack(maxObjects+1);
    searchStack.push(weightEntry);
    
    //---

    while(!searchStack.empty())
    {   
        vector<Weight> allowableWeights,curWeights;
        
        unsigned insertOffset;

        //Check the top of the stack
        //--
        insertOffset=searchStack.top().offset;

    
        //Don't proceed if we have used all the child weights   
        if(insertOffset >=searchStack.top().allowableWeights.size())
        {
            //We have run out of allowable weights to use
            searchStack.pop();
            continue;
        }

        //Tell current stack to move to next child
        searchStack.top().offset++;

        //Copy the top of the stack
        curWeights=searchStack.top().curWeights;
        allowableWeights=searchStack.top().allowableWeights;
        
        //Insert a weight into the sack
        float remainingWeight,currentWeight;
        curWeights.push_back(allowableWeights[insertOffset]);
        currentWeight=searchStack.top().cumulativeWeight +allowableWeights[insertOffset].mass;

        for(unsigned int ui=0;ui<targetWeight.size();ui++)
        {
            if( currentWeight <= targetWeight[ui]+tolerance &&  
                    currentWeight >=targetWeight[ui]-tolerance)
            {
                //TODO: Use solution tree, to more efficiently store possible solutions?
                solutions.push_back(curWeights);
            }
        }
            

        //Erase disallowed weights.
        //the weights are sorted in increasing size, so just
        // cut once we hit the right threshold
        //--------
        remainingWeight=upperWeight-currentWeight;
        if(remainingWeight > 0)
        {
            for(size_t ui=0;ui<allowableWeights.size();ui++)
            {
                if(allowableWeights[ui].mass > remainingWeight)
                {
                    allowableWeights.resize(ui+1);
                    break;
                }
            }
        }   
        else 
            allowableWeights.clear();
        //--------

        //Add the next weight to the stack, if we have not run out of number of objects
        if(allowableWeights.size() && curWeights.size() != maxObjects)
        {
            weightEntry.curWeights.swap(curWeights);
            weightEntry.allowableWeights.swap(allowableWeights);
            weightEntry.cumulativeWeight=currentWeight;
            //Search in decreasing weight only, to not generate duplicates
            weightEntry.offset=insertOffset;
            searchStack.push(weightEntry);
        }


    }

}

#ifdef DEBUG
bool MassTool::runUnitTests()
{
    vector<Weight> weights;
    weights.push_back(1.0);
    weights.push_back(3.0);
    weights.push_back(10.0);


    //Test the knapsack algorithm. Remember some input
    // weights may be removed during this calculation
    vector<vector<Weight> > solutions;
    bruteKnapsack(weights,4.0f,0.01f,4,solutions);

    TEST(solutions.size() ==2,"Brute-force knapsack test");

    solutions.clear();
    weights.clear();

    //set up the problem again
    weights.push_back(1.0);
    weights.push_back(3.0);
    weights.push_back(10.0);

    vector<float> targetWeights;
    targetWeights.push_back(4.0f);
    targetWeights.push_back(10.0f);
    bruteKnapsack(weights,targetWeights,
                0.01f,4,solutions);

    TEST(solutions.size() ==4,"Brute-force knapsack test with multiple targets");
    return true;    
}

#endif