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
>>> 

A polynomial time exact algorithm for the Subset sum problem

# ssp.py - Polynomial time algorithm for Subset sum problem.
#
# Copyright 2020 Andrea Bianchini 
#
# Theory of the implemented algorithm can be found at:
#                                   http://es-andreabianchini.it/ssp_eng.pdf
#
#  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/>.
#

from random import randrange
import math


def computeBase(c12):
    i=0;
    while(c12[i]==0 and i<n-1):
        i=i+1
    # i is the first one
    i=i+1
    # skipped
    while(i<n-1 and c12[i]==0):
    	i=i+1
    # i is the second one
    while(i<n-1 and c12[i]==1):
    	i=i+1

    if (i>n-1):
    	return n-1

    return i

def findMax(topvalue, currvalue, value, start):
    lsb1=0
    lsb2=0
    lsb22=0
    lsbmax=0
    lsbmin=0
    i=0
    basemax=0
    base=0
    m=0
    firstbit=0
    k=0
    amax=0
    kk=0
    oldamax=0
    solfound=0
    global ro
    global teta
    global alfa
    global niter
    global niter1
    global w
    global sol
    global c1
    global bestsolk
    
    if (ro<0.0):
            ro=value*2.0/(w[start]+w[0])
    #       printf("ro=%.3f\n",ro)


    if (teta<0.0):
            i=1
            teta=math.log(n)/math.log(10.0)
            log10N=teta
            while(log10N>i):
                    teta*=math.log(n/(math.pow(10.0,log10N-(i))))/math.log(10.0)
                    i=i+1

            teta/=2.0
            teta=teta+value*2.0/(n*(w[n-1]+w[0]))
            teta=((teta*10.0))/10.0
    #   printf("teta=%.3f\n",teta)


    niter1=niter1+1

    kk=0
    for i in range(start):
            kk=kk+w[i]

    if (currvalue+kk<=bestsolk):
            return

    alfa=alfa+1

    i=start
    k=0;
    lsb1=-1
    firstbit=-1
    while(k<value and i>=0):
            niter=niter+1
            if (k+w[i]<=value):
                    k=k+w[i]
                    lsb1=i
                    if (firstbit==-1):
                            firstbit=i
            else:
                    if (lsb1>-1):
                            break
            i=i-1
            
    lsbmin=0
    lsbmax=lsb1-1
    while(lsbmin+1<lsbmax):
            niter=niter+1
            m=(lsbmin+lsbmax)//2
            if (k+w[m]>value):
                    lsbmax=m
            else:
                    lsbmin=m

    lsb2=lsbmax
    if (lsb2>-1 and k+w[lsb2]<=value and  lsb2<lsb1):
            k=k+w[lsb2]
    else:
            lsb2=lsbmin
            if (lsb2>-1 and k+w[lsb2]<=value and lsb2<lsb1):
                    k=k+w[lsb2]
            else:
                    lsb2=-1

    #    printf("value-k=%I64d\n",value-k);
    amax=k
    oldamax=amax

    basemax=firstbit+1
    base=firstbit+1

    while(base>1):
    #       printf("base=%d LSB2=%d\n",base,lsb2)
    #       getch()
            if (lsb2>=0):
                    k=k-w[lsb2]
            i=base-1
            if (i>-1):
                    k=k-w[i]
            i=lsb1-1
            lsb1=-1
            while(k<value and i>=0):
                    niter=niter+1
                    if (k+w[i]<=value):
                            k=k+w[i]
                            lsb1=i
                    else:
                            if (lsb1>-1):
                                    break
                    i=i-1

            lsbmin=0
            lsbmax=lsb1-1
            while(lsbmin+1<lsbmax):
                    niter=niter+1
                    m=(lsbmin+lsbmax)//2
                    if (k+w[m]>value):
                            lsbmax=m
                    else:
                            lsbmin=m

            lsb2=lsbmax
            if (lsb2>=0 and k+w[lsb2]<=value and lsb2<lsb1):
                    k=k+w[lsb2]
            else:
                    lsb2=lsbmin
                    if (lsb2>=0 and k+w[lsb2]<=value and lsb2<lsb1):
                            k=k+w[lsb2]
                    else:
                            lsb2=-1

            nitems=base-1-lsb1
            if lsb2>-1:
                    nitems=nitems+1
            if (k>amax and lsb2>=0 and lsb2<lsb1 and nitems<=ro):
                    oldamax=amax
                    amax=k
                    basemax=base-1

            base=base-1

    fbit=-1
    sbit=-1
    i=basemax-1
    k=0
    lsb1=-1
    while(k<value and i>=0):
            niter=niter+1
            if (k+w[i]<=value):
                    k=k+w[i]
                    lsb1=i
                    c1[i]=1
                    if (fbit==-1):
                            fbit=i
                    if (fbit>-1 and sbit==-1):
                            sbit=i
            else:
                    if (lsb1>-1):
                            break
            i=i-1
            
    lsbmin=0
    lsbmax=lsb1-1
    while(lsbmin+1<lsbmax):
            niter=niter+1
            m=(lsbmin+lsbmax)//2
            if (k+w[m]>value):
                    lsbmax=m
            else:
                    lsbmin=m

    lsb2=lsbmax
    if (lsb2>-1 and k+w[lsb2]<=value and lsb2<lsb1):
            k=k+w[lsb2]
            c1[lsb2]=1
    else:
            lsb2=lsbmin
            if (lsb2>-1 and k+w[lsb2]<=value and lsb2<lsb1):
                    k=k+w[lsb2]
                    c1[lsb2]=1
            else:
                    lsb2=-1

    if (sbit<0):
            sbit=lsb2

    k1=currvalue+k

    if (k1>bestsolk):
            print("basemax="+str(basemax)+" sol="+str(k1))
            bestsolk=k1
            for zz in range(n):
                    sol[zz]=c1[zz]

    if (k==value):
            solfound=1
            alfa=alfa-1
            return

    superbase=computeBase(c1)
    lsb=0

    amenaprimo=0;
    bprimo=0;
    b=-1;
    cbase = superbase;
    a=k;
    k=value-k;
    cnt=0;

    if (lsb1==-1):
            alfa=alfa-1
            return

    lsb=0
    while(c1[lsb]==0 and lsb<n):
            lsb=lsb+1
    if (lsb>=n):
            lsb=basemax+1
    lsb=lsb-1

    alfa1=1


    if (superbase<=start+1):
            while(solfound==0 and lsb>-1 and lsb<superbase-3 and teta*bprimo<=a-amenaprimo and c1[fbit]==1 and c1[sbit]==1):
                    cnt=cnt+1
                    if (lsb>-1):
                            findMax(topvalue,topvalue-k,k,lsb)

                    for i1 in range(lsb+1):
                            c1[i1]=0

                    while(c1[lsb]==0 and lsb<superbase):
                            lsb=lsb+1
                    if (lsb>=superbase):
                            lsb=superbase-1
                    c1[lsb]=0
                    k=k+w[lsb]
                    amenaprimo=amenaprimo+w[lsb]
                    bprimo=value-(a-amenaprimo)
                    lsb-=1

    alfa=alfa-1
    return

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)]
    ro=-1.0
    teta=-1.0
    niter=0
    niter1=0
    bestsolk=0
    alfa=0

    print("n="+str(n))
    print("c="+str(c))
    print("wmin="+str(wmin))
    print("wmax="+str(wmax))
    print("w = ",w)
    print()

    findMax(c,0,c,n-1)
    solvalue=sum([sol[i]*w[i] for i in range(n)])
    print()
    print("Solution="+str(solvalue))
    print("Sol. Items = ",[w[i] for i in range(n) if sol[i]==1])
    
except Exception:
    print()
    print("ssp by Andrea Bianchini.")
    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()