An Exact Algorithm of Polynomial Complexity for the Subset Sum Problem

# ssp.py - An Exact Algorithm of Polynomial Complexity for the Subset Sum Problem
#
# Andrea Bianchini, 2024, Common Creative CC BY-NC 4.0. 
#
# see also https://es-andreabianchini.it/ssp_eng.pdf
#

from random import randrange
import math

try:
    inp = open("input.txt","rt")
    n=int(inp.readline())
    c=int(inp.readline())
    wmin=int(inp.readline())
    wmax=int(inp.readline())
    inp.close()
    w=sorted([int(wmin+randrange(wmax-wmin)) for i in range(n)])
    c1=[0 for i in range(n)]
    sol=[0 for i in range(n)]
    bestsolu=[0 for i in range(n)]    
    niter=0
    bestsol=0
    bestsolk=0

except Exception:
    print(Exception)
    print()
    print("ssp, Andrea Bianchini, 2024, Common Creative CC BY-NC 4.0.")
    print("USAGE:")
    print("ssp")
    print("ssp reads input data from input.txt")
    print("WHERE input.txt is in the form:")
    print("n     # number of items")
    print("c     # capacity of the knapsack")
    print("wmin  # min item's weight")
    print("wmax  # max item's weight")
    print()
    
def solve():
    global n
    global c
    global w
    global bestsolu
    global sol
    global bestsol
        
        
    bestsol=0
    solvalue=0
    divi=1
    T=int(pow(2,n-1))
    while(T>1 and bestsol<c):
        t=0
        divi=1
        while(divi<int(pow(2,n)) and bestsol<c):

            #print(bin(T),bin(t),bestsol)


            sol=list(map(int,bin(t+int(T/divi))[2:].zfill(n)[::-1]))
            gt=sum([sol[i]*w[i] for i in range(n)])
            if (gt<=c):
                t=t+int(T/divi)

            divi*=2
            
            solvalue=gt
            if (bestsol<solvalue and solvalue<=c):
                bestsolu=list(map(int,bin(t)[2:].zfill(n)[::-1]))
                bestsol=solvalue

        T=int(T/2)
                   
    return

try:
    while(True):
        w=sorted([int(wmin+randrange(wmax-wmin)) for i in range(n)])
        
        print("n="+str(n))
        print("c="+str(c))
        print("wmin="+str(wmin))
        print("wmax="+str(wmax))
        print("w = ",w)
        print()

        bestsol=0
        
        solve()
    
        solvalue=sum([bestsolu[i]*w[i] for i in range(n)])

        print()
        print("Solution="+str(solvalue))
        print("Sol. Vector = ",bestsolu)
        print("Sol. Items = ",[w[i] for i in range(n) if bestsolu[i]==1])
        print("Hit return...")
        input()
except Exception:
    print(Exception.message)

input.txt file example:

50
10000
200
4000

Execution example:

n=50
c=10000
wmin=200
wmax=4000
w =  [257, 537, 580, 580, 664, 671, 688, 755, 844, 939, 953, 960, 1014, 1095, 1250, 1359, 1417, 1463, 1513, 1527, 1575, 1629, 1737, 1738, 1817, 1841, 1935, 1942, 2115, 2240, 2318, 2330, 2565, 2606, 2681, 2692, 2788, 3075, 3112, 3196, 3201, 3230, 3257, 3368, 3400, 3440, 3508, 3512, 3782, 3938]


Solution=10000
Sol. Vector =  [1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
Sol. Items =  [257, 580, 580, 664, 671, 688, 755, 844, 939, 953, 960, 1014, 1095]
Hit return...

n=50
c=10000
wmin=200
wmax=4000
w =  [293, 419, 439, 576, 635, 672, 978, 1031, 1063, 1120, 1229, 1285, 1288, 1446, 1660, 1800, 1978, 2014, 2033, 2034, 2040, 2148, 2185, 2193, 2201, 2479, 2506, 2521, 2540, 2553, 2585, 2673, 2799, 2976, 3051, 3089, 3287, 3371, 3379, 3390, 3438, 3581, 3613, 3661, 3696, 3753, 3763, 3766, 3834, 3836]


Solution=10000
Sol. Vector =  [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
Sol. Items =  [293, 2201, 2479, 2506, 2521]
Hit return...

https://es-andreabianchini.it/ssp_eng.pdf

Andrea Bianchini 2024.

Andrew’s genetic algorithm

In questo articolo vi presento un possibile algoritmo genetico da me implementato in C++ per questioni di velocità di calcolo (e non solo…).

L’algoritmo è molto semplice e compatto, ma vi assicuro che funziona benissimo! Come esempio di problema da risolvere ho adottato il subset sum problem.

//==========================================================================
// Name        : aga.cpp
// Author      : Andrea Bianchini (2021)
// Version     : 1.0.0
// Copyright   : SOME RIGHT RESERVED CREATIVE COMMONS BY 4.0
// License     : https://creativecommons.org/licenses/by/4.0/
// Description : AGA - Andrew's Genetic Algorithm
// Use         : AGA p = AGA(PSIZE,ISIZE,RMIN,RMAX,IMAX,PCROSS,PMUT);
//             : p.evolve();
// Where       : PSIZE = Population size (number of individuals)
//             : ISIZE = Individual size (number of variables)
//             : RMIN  = Individual variables min value
//             : RMAX  = Individual variables max value
//             : IMAX  = Max number of iteration
//             : PCROSS= Crossover probability
//             : PMUT  = Mutation probability
// Example     : AGA p = AGA(500,20,0,1,10,0.5,0.3);
//             : p.evolve();
//             :
//             : Personalize your objective into evalIndividual().
//             : This file is personalized to solve subset sum problem.
//==========================================================================
#include <stdlib.h>
#include <limits.h>
#include <stdio.h>
#include <math.h>
#include <time.h>


class AGA
{
protected:
	int POP_SIZE = 500;
	int IND_SIZE = 20;
	int RANGE_MIN = 0;
	int RANGE_MAX = 1;
	int MAXITER = 10;
	float PCROSS = 0.5;
	float PMUT = 0.3;
	int **pop;
	int first=1;
	int solfound;
	int niter;

public:
	int *w;
	int *bestIndividual;
	int c = 400;

	void initSSP()
	{
		w=new int [IND_SIZE];
		for(int i=0;i<IND_SIZE;i++)
		w[i]=c*0.1+(int)(rand()%(int)(c*0.3-c*0.1+1));

	}
	AGA(int pops, int inds, int rmin, int rmax, int maxiter, float pcross, float pmute)
		{
			solfound=0;
			niter=0;
			first=1;
			POP_SIZE=pops;
			IND_SIZE=inds;
			RANGE_MIN=rmin;
			RANGE_MAX=rmax;
			MAXITER = maxiter;
			PCROSS = pcross;
			PMUT = pmute;
			pop = new int* [POP_SIZE];
			int i;
			for(i=0;i<POP_SIZE;i++)
			{
				pop[i] = new int [IND_SIZE+1]; // +1 for fit
				pop[i][IND_SIZE]=0;
			}
			bestIndividual=new int [IND_SIZE+1];
			bestIndividual[IND_SIZE]=0;
			initSSP();
		}

	~AGA()
	{
		for(int i=0;i<POP_SIZE;i++)
			delete[] pop[i];
		delete[] pop;
		delete[] bestIndividual;
		delete[] w;
	}

	int * individual(int i)
	{
		return pop[i];
	}

	void populate()
	{
		int ind;
		for(ind=0;ind<POP_SIZE-1;ind++)
			for(int k=ind+1;k<POP_SIZE;k++)
				if (abs(individual(ind)[IND_SIZE]-c)>abs(individual(k)[IND_SIZE]-c))
				{
					int *p=individual(ind);
					pop[ind]=pop[k];
					pop[k]=p;
				}
		//..........
		if (first==1)
		{
			for(ind=0;ind<POP_SIZE;ind++)
			{
				int *in;
				in = individual(ind);
				for(int j=0;j<IND_SIZE;j++)
				{
					in[j] = rndInRange();
				}
				in[IND_SIZE]=0; // fit
			}
			first=0;
		}
		else
		{
			int p2=POP_SIZE*PCROSS+1;
			int p1=0;
			while(p2<POP_SIZE && p1<POP_SIZE*PCROSS+1)
			{
				Cross(p1,p2);
				if (p2<POP_SIZE*PMUT)
					Mutate(p2);
				p1++;
				p2++;
			}
			while (p2<POP_SIZE)
			{
				int *in;
				in = individual(p2);
				for(int j=0;j<IND_SIZE;j++)
				{
					in[j] = rndInRange();
				}
				in[IND_SIZE]=0; // fit
				p2++;
			}
		}
	}

	void Cross(int p1, int p2)
	{
		int *c1,*c2,*c3;

		c1=individual(p1);
		c2=individual(p1+1);
		c3=individual(p2);

		for(int i=0;i<IND_SIZE/2;i++)
		{
			c3[i]=c1[i];
		}
		for(int i=IND_SIZE/2;i<IND_SIZE;i++)
		{
			c3[i]=c2[i];
		}
		c3[IND_SIZE]=0;
	}

	void Mutate(int p2)
	{
		int c = rand() % IND_SIZE;
		int *d=individual(p2);
		d[c] = rndInRange();
	}

	int rndInRange()
	{
		int a = RANGE_MIN+(int)(rand()%(RANGE_MAX-RANGE_MIN+1));
		if (a>RANGE_MAX)
			a=RANGE_MAX;
		//printf("%d\n\r",a);
		return a;
	}

	void evalPopulation()
	{
		for(int i=0;i<POP_SIZE;i++)
		{
			evalIndividual(i);
		}

	}

	void evalIndividual(int i)
	{
		int z=0;
		int *p = individual(i);
		for(int j=0;j<IND_SIZE;j++)
			z+=w[j]*p[j];
		//if (z<c)
			p[IND_SIZE]=z;
		//else
			//p[IND_SIZE]=0;
		if (z<=c && z>bestIndividual[IND_SIZE])
		{
			for(int j=0;j<IND_SIZE;j++)
				bestIndividual[j]=p[j];
			bestIndividual[IND_SIZE]=z;
			printf("Best Fit=%d\r",z);
			if (z==c)
			{
				solfound=1;
				printf("Optimal Solution z=%d found!\r",z);
			}
		}
	
	}
	
	void evolve()
	{
		niter=0;
		while(niter<MAXITER && solfound==0)
		{
			printf("--------------------\rGenerazione n.%d\r",niter);
			populate();
			evalPopulation();
			int ind;
			int tot=0;
			for(ind=0;ind<POP_SIZE-1;ind++)
			{
				tot+=individual(ind)[IND_SIZE];
				for(int k=ind+1;k<POP_SIZE;k++)
					if (individual(ind)[IND_SIZE]<individual(k)[IND_SIZE])
					{
						int *p=individual(ind);
						pop[ind]=pop[k];
						pop[k]=p;
					}
			}
			tot+=individual(POP_SIZE-1)[IND_SIZE];
			tot/=POP_SIZE;
			printf("MAX=%d   AVG=%d   MIN=%d\r",pop[0][IND_SIZE],tot,pop[POP_SIZE-1][IND_SIZE]);
			niter++;
		}
	}
};

int main() {
	int POP_SIZE=500;
	int IND_SIZE=20;
	int RANGE_MIN=0;
	int RANGE_MAX=1;
	int MAXITER = 10;
	float PCROSS = 0.5;
	float PMUT = 0.3;

	AGA p = AGA(POP_SIZE,IND_SIZE,RANGE_MIN,RANGE_MAX,MAXITER,PCROSS,PMUT);
	p.c=600;
	p.initSSP();
	clock_t tStart = clock();
	p.evolve();
	printf("\rTempo di esecuzione: %.4fs\r", (double)(clock() - tStart)/CLOCKS_PER_SEC);

	printf("\r---------------\rSoluzione: %d\r",p.bestIndividual[IND_SIZE]);
	printf("Items:\r[");
	for(int i=0;i<IND_SIZE;i++)
		printf("%d ",p.w[i]);
	printf("]\r");
	printf("Items soluzione:\r[");
	for(int i=0;i<IND_SIZE;i++)
		printf("%d ",p.w[i]*p.bestIndividual[i]);
	printf("]\r");
	return 0;
}

Esempio di output:

--------------------
Generazione n.0
Best Fit=524
Best Fit=548
Best Fit=579
MAX=1833   AVG=382   MIN=339
--------------------
Generazione n.1
Best Fit=593
MAX=1689   AVG=333   MIN=283
--------------------
Generazione n.2
Best Fit=598
MAX=1189   AVG=351   MIN=283
--------------------
Generazione n.3
Best Fit=600
Optimal Solution z=600 found!
MAX=1189   AVG=406   MIN=203

Tempo di esecuzione: 0.0267s

---------------
Soluzione: 600
Items:
[92 123 109 82 174 110 148 147 96 136 169 67 114 129 100 124 122 101 62 81 ]
Items soluzione:
[92 0 0 0 0 0 0 0 0 136 0 67 0 0 100 124 0 0 0 81 ]

Coppie di numeri a differenza costante

# kdiff.py
#
# 
# Questo programma risolve il seguente problema :
# Dato un insieme s di N numeri interi trovare
# tutte le coppie di interi la cui differenza è pari a K
#
# by Andrea Bianchini (2021)

from random import randint

N = 50
MIN = 1
MAX = 1000
K = randint(MIN,MAX/2)

s = [randint(MIN,MAX) for _ in range(N)]


def f(x,y):
    global K
    if abs(x-y)==K:
        return True
    else:
        return False
    
res =[]

for x in range(N-1):
    for y in range(x+1,N):
        if f(s[x],s[y]):
            res.append((s[x],s[y]))

print("s =",s)
print()
print("K = %d" %K)
print("coppie =",res)

Esempio 1:

s = [240, 662, 961, 108, 802, 232, 954, 99, 542, 775, 668, 988, 637, 96, 135, 549, 111, 678, 86, 596, 912, 13, 191, 516, 118, 88, 837, 632, 870, 843, 542, 569, 819, 579, 411, 154, 743, 89, 95, 686, 604, 503, 511, 859, 305, 671, 161, 510, 164, 906]

K = 10
coppie = [(108, 118), (99, 89), (668, 678), (96, 86), (569, 579), (154, 164)]
>>> 

Esempio 2:

s = [327, 847, 660, 722, 431, 287, 79, 542, 354, 514, 703, 522, 764, 75, 904, 789, 103, 801, 929, 884, 964, 818, 992, 458, 369, 900, 239, 671, 120, 904, 421, 342, 403, 314, 495, 794, 119, 140, 810, 237, 333, 639, 676, 705, 134, 802, 528, 20, 452, 2]

K = 105
coppie = [(239, 134), (342, 237), (810, 705)]
>>> 

Algoritmi genetici

# DEAP_generico.py
#
# esempio di utilizzo della libreria per la computazione genetica evolutiva
# dei problemi denominata DEAP.
# in questo caso ho utilizzato lo scheletro della risoluzione di un
# problema generico per risolvere il problema di ottimizzazione
# noto come Knapsack Problem.
# documentazione su : https://deap.readthedocs.io/en/master/index.html
#
# by Andrea Bianchini (2021)



import random

from deap import base
from deap import creator
from deap import tools

IND_SIZE = 50
C = 10000
MIN_WEIGHT = int(C*0.01)
MAX_WEIGHT = int(C*0.07)
MAX=0
sol=[]
items=[random.randint(MIN_WEIGHT,MAX_WEIGHT) for _ in range(IND_SIZE)]

creator.create("FitnessMax", base.Fitness, weights=(1.0,))
creator.create("Individual", list, fitness=creator.FitnessMax)

toolbox = base.Toolbox()
# Attribute generator 
toolbox.register("attr_bool", random.randint, 0, 1)
# Structure initializers
toolbox.register("individual", tools.initRepeat, creator.Individual, 
    toolbox.attr_bool, IND_SIZE)
toolbox.register("population", tools.initRepeat, list, toolbox.individual)

def evalOneMax(individual):
    global sol
    global MAX
    e1 = sum([items[i]*individual[i] for i in range(len(individual))])
    if e1>MAX and e1<=C:
        MAX = e1
        sol=individual
    return e1,

toolbox.register("evaluate", evalOneMax)
toolbox.register("mate", tools.cxTwoPoint)
toolbox.register("mutate", tools.mutFlipBit, indpb=0.05)
toolbox.register("select", tools.selTournament, tournsize=3)

summa=0
def main():
    global summa
    pop = toolbox.population(n=300)

    # Evaluate the entire population
    fitnesses = list(map(toolbox.evaluate, pop))
    for ind, fit in zip(pop, fitnesses):
        ind.fitness.values = fit

    # CXPB  is the probability with which two individuals
    #       are crossed
    #
    # MUTPB is the probability for mutating an individual
    CXPB, MUTPB = 0.5, 0.2

    # Extracting all the fitnesses of 
    fits = [ind.fitness.values[0] for ind in pop]

    # Variable keeping track of the number of generations
    g = 0
    
    # Begin the evolution
    
    while max(fits) < C and g < 1000:
        # A new generation
        g = g + 1
        print("-- Generation %i --" % g)

        # Select the next generation individuals
        offspring = toolbox.select(pop, len(pop))
        # Clone the selected individuals
        offspring = list(map(toolbox.clone, offspring))

        # Apply crossover and mutation on the offspring
        for child1, child2 in zip(offspring[::2], offspring[1::2]):
            if random.random() < CXPB:
                toolbox.mate(child1, child2)
                del child1.fitness.values
                del child2.fitness.values

        for mutant in offspring:
            if random.random() < MUTPB:
                toolbox.mutate(mutant)
                del mutant.fitness.values

        # Evaluate the individuals with an invalid fitness
        invalid_ind = [ind for ind in offspring if not ind.fitness.valid]
        fitnesses = map(toolbox.evaluate, invalid_ind)
        for ind, fit in zip(invalid_ind, fitnesses):
            ind.fitness.values = fit

        pop[:] = offspring

        # Gather all the fitnesses in one list and print the stats
        fits = [ind.fitness.values[0] for ind in pop]
        
        length = len(pop)
        mean = sum(fits) / length
        sum2 = sum(x*x for x in fits)
        std = abs(sum2 / length - mean**2)**0.5
        
        print("  Min %s" % min(fits))
        print("  Max %s" % max(fits))
        print("  Avg %s" % mean)
        print("  Std %s" % std)

        if max(fits) <=C and summa<max(fits):
            summa = max(fits)

main()
print()
print("Soluzione ottima = %d" %MAX)
print("Lista items completa :")
print(items)
print("Lista items soluzione :")
osol = [items[i]*sol[i] for i in range(IND_SIZE)]
print(osol)

Esempio:

Soluzione ottima = 10000
Lista items completa :
[481, 167, 361, 111, 670, 602, 119, 331, 654, 155, 400, 227, 651, 531, 413, 134, 409, 140, 488, 457, 495, 628, 153, 632, 327, 159, 593, 633, 682, 392, 368, 526, 599, 478, 408, 315, 466, 582, 302, 172, 427, 173, 551, 673, 272, 158, 502, 269, 685, 588]
Lista items soluzione :
[0, 0, 361, 0, 670, 0, 0, 331, 654, 0, 0, 0, 651, 0, 0, 0, 409, 140, 488, 0, 0, 628, 0, 632, 327, 0, 593, 0, 682, 0, 368, 0, 0, 478, 0, 315, 0, 582, 0, 172, 427, 0, 551, 0, 272, 0, 0, 269, 0, 0]
>>> 

Programmazione Lineare

# pulp1.py
#
# Risoluzione di un problema di programmazione lineare
# tramite la libreria PuLP.
#
# definizione problema by https://www3.diism.unisi.it/~agnetis/esesvolti.pdf
# soluzione in Python tramite utilizzo libreria PuLP, by Andrea Bianchini 2021
#
# Problema:
# Un lanificio produce filato di tipo standard e di tipo speciale
# utilizzando 3 diverse macchine, le cui produzioni orarie sono le seguenti:
# macchina A: 3 matasse standard e 1 speciale
# macchina B: 2 matasse standard e 2 speciali
# macchina C: 2 matasse standard e 1 speciale
# Il mercato richiede almeno 60 matasse standard e 40 di tipo speciale al giorno. I costi
# orari delle due macchine sono: 90 euro per la A, 80 euro per B, 60 euro per C.
# Scrivere un modello di programmazione lineare per determinare la produzione giornaliera
# di costo minimo. (Non occorre imporre il vincolo che le ore giornaliere non superino 24)
#


from pulp import *

a = pulp.LpVariable("a", lowBound=0)
b = pulp.LpVariable("b", lowBound=0)
c = pulp.LpVariable("c", lowBound=0)

problem = pulp.LpProblem("Un semplice problema di min", LpMinimize)

problem += 90*a + 80*b + 60*c, "The objective function"
problem += 3*a + 2*b + 2*c >= 60, "1st constraint"
problem += a + 2*b + c >= 40, "2nd constraint"
problem += a >= 0, "3rd constraint"
problem.solve()

print("Risultati della ottimizzazione:")
for variable in problem.variables():
    print(variable.name + "=" + str(variable.varValue))
print("Costo minimo netto totale: %.1f" %value(problem.objective))

Esempio:

Risultati della ottimizzazione:
a=0.0
b=10.0
c=20.0
Costo minimo netto totale: 2000.0
>>>