Tävlingsprogrammering/Uppgifter/Isomerer

Isomerer

En molekyl av en viss sorts alkoholer kan beskrivas på följande sätt:

Den innehåller totalt N kol-atomer.
Exakt en av dessa är bunden till en s.k. hydroxylgrupp (–OH). Denna kol-atom kan binda till högst tre andra kol-atomer.
Varje övrig kol-atom kan binda till högst fyra andra kol-atomer.
Alla kol-atomer måste vara sammanlänkade.
Inga cykler får förekomma.

Alla isomerer för N=1, N=2, N=3 och N=5. För N=4 finns fyra isomerer. Varje kula motsvarar en kol-atom och varje streck en bindning.

De enklaste av sådana molekyler är metanol (N = 1) och etanol (N = 2; “vanlig” alkohol). Från och med N = 3 finns det flera olika molekyler som uppfyller kraven ovan men som har olika struktur (se figur ovan). Du ska skriva ett program som frågar efter antalet kolatomer N , där 3 ≤ N ≤ 11, och beräknar hur många olika godkända molekyler, s.k. isomerer, det finns.

Två molekyler räknas som olika om det inte går att numrera kolatomerna i vardera molekylen på ett sådant sätt att samma par av kol-atomer är sammanbundna i båda molekylerna och att samma kol-atom är bunden till hydroxylgruppen.

Denna numrering av kol-atomerna avslöjar att detta är två avbildningar av samma isomer.

Körningsexempel:

Antal kolatomer ? 6  
Antal isomerer: 17

Lösning[redigera]

Detta är ett klassiskt problem som löstes redan år 1874 ^[1]. Att fråga efter antalet alkoholer (med OH-grupp) istället för antalet alkaner (utan OH-grupp) kan verka som en extra svårighet, men det är faktiskt precis tvärtom, det gör problemet betydligt lättare. Anledningen är att man har en fast punkt (kolatomen där OH-gruppen sitter) att utgå ifrån när man t.ex. ska avgöra om två molekyler är olika. I grafteoretiska termer handlar det om antalet "rooted trees" med ett visst antal noder istället för antalet "unrooted trees".

Det går att hitta en hyfsat enkel rekursionsformel som ger svaret. Den är dock inte helt lätt att komma fram till, så gränsen på uppgiften är vald så att man även ska kunna "rita upp" varje tänkbar molekyl (bildligt talat) och ta bort dubletter. Båda lösningarna beskrivs här.

Rekursionsformel[redigera]

Rekursionsformeln lyder så här ^[2]:

$R(0)=1\,$

$R(1)=1\,$

$R(n+1)=\sum _{a<b<c}^{a+b+c=n}R(a)R(b)R(c)+\sum _{a\neq c}^{2a+c=n}{{R(a)+1} \choose 2}R(c)+\sum ^{3a=n}{{R(a)+2} \choose 3}$

där summorna ska gå över de kombinationer av icke-negativa heltal a, b och c som uppfyller kraven (sista termen förekommer alltså enbart när n är delbart med 3). Det är ganska lätt att förstå vad de tre termerna betyder. För att räkna ut R(n+1) utgår man från "roten" i trädet (den kolatom som har OH-gruppen). De n övriga kolatomerna kan vara fördelade på tre grenar med olika antal i varje (första termen), med två lika och den tredje olika (andra termen) eller med alla tre lika antal kolatomer (tredje termen). Grenar kan förstås även vara tomma.

Rekursionsformeln går att implementera ganska direkt. Om man dessutom sparar de beräknade värdena på R(n) så blir programmet supersnabbt (dynamisk programmering).

#include <stdio.h>

int over(int n,int k) {
  int s=1,i;
  for(i=0;i<k;i++) s*=n-i;
  for(i=0;i<k;i++) s/=i+1;
  return s;
}

int main() {
  int n, R[100],i,j,k,N;
  printf("Antal kolatomer ?");
  scanf("%d", &N);
  R[0]=R[1]=1;
  for(n=2;n<=N;n++) {
    R[n]=0;
    for(i=0;3*i+3<=n-1;i++) for(j=i+1;i+j+j+1<=n-1;j++) {
        k=n-1-i-j;
        R[n]+=R[i]*R[j]*R[k];
      }
    for(i=0;2*i<=n-1;i++) {
      j=n-1-2*i;
      if(i!=j) R[n]+=over(R[i]+1,2)*R[j];
    }
    if((n-1)%3==0) R[n]+=over(R[(n-1)/3]+2,3);
  }
  printf("Antal isomerer: %d\n", R[N]);
  return 0;
}

Utan rekursionsformeln[redigera]

Man behöver så klart inte kunna rekursionsformeln sedan tidigare för att räkna isomererna och det behöver inte nödvändigtvis vara krångligare. Anta att vi har ett träd T med n noder, rot r och som har högst 3 delträd T_1, T_2, T_3 som alla tre också har samma struktur. Då är

$|T_{1}|+|T_{2}|+|T_{3}|=n-1$

Vi antar att

$|T_{1}|=a,|T_{2}|=b,|T_{3}|=c,a\leq b\leq c$

Då måste a+b+c = n-1 och $\sigma =\langle a,b,c\rangle$ är en ickeavtagande partition av n-1. Kan vi beräkna antalet träd om vi vet σ? Ja, man behöver bara känna till R(a), R(b), R(c) och akta sig för dubbelräkning.

Här är programmet i python: Det genererar först alla partitioner (a, b, c) och räknar antalet träd för varje partition.

def partition(n, m):
    if m == 1:  yield [n]
    elif m > 1:
        for i in range(n//m + 1):
            for P in partition (n - m*i, m-1):
                yield [i] + [x+i for x in P]

def R(n):
    if n <= 2: return 1
    summa = 0
    for (a, b, c) in partition(n-1, 3):
        A, B, C = R(a), R(b), R(c) # för att inte räkna dem igen
        if a == b and c != a:
            summa += A*(A-1)*C//2 # fallet alla olika, (A choose 2) * C
            summa += A*C # fallet a-trädet isomorft med b-trädet
        elif b == c and c != a:
            summa += B*(B-1)*A//2 # alla olika, (B choose 2)*A
            summa += B*A # B, C -träd isomorfa
        elif a == c:
            summa += A*(A-1)*(A-2)//6 # alla olika, (A choose 3)
            summa += A*(A-1) # 2 isomorfa, (A choose 2) * 2
            summa += A # alla isomorfa
        else:
            summa += A*B*C # alla måste vara olika
    return summa
n = input("n? ")
print("R(n) = ", R(n))

"Rita upp" molekylerna[redigera]

Det finns många sätt att göra detta. Programmet blir betydligt krångligare än med rekursionsformeln, men det kan kännas "säkrare" att skriva. Följande lösning använder sig av en representation av varje moleyl som en array av heltal, där varje heltal k anger hur många grenar kolatomen har (borträknat den man kom ifrån) och de följande talen beskriver de k grenarna enligt samma princip. Första talet i arrayen betecknar kolatomen som binder OH-gruppen och man tänker sig att man "kommer från" denna. För den som är kemiskt bildad kan följande exempel hjälpa:

Molekylens namn	Representation
Metanol	0
Etanol	1,0
1-propanol	1,1,0
2-propanol	2,0,0
1-butanol	1,1,1,0
2-butanol	2,1,0,0 (eller 2,0,1,0)

För varje N går man igenom alla molekyler med N-1 kolatomer och lägger till en atom på alla möjliga ställen samt kollar om den skapade molekylen är en dublett, detta kräver lite eftertanke.

#include <stdio.h> 
#include <stdlib.h>

int blen(int *p) {
  int i,t=1;
  for(i=0;i<p[0];i++) t+=blen(p+t);
  return t;
}

int equal(int *p, int *r) {
  int p1,p2,p3,r1,r2,r3;
  if(p[0]!=r[0]) return 0;
  if(p[0]==0) return 1;
  if(p[0]==1) return equal(p+1,r+1);
  if(p[0]==2) {
    p1=blen(p+1);
    r1=blen(r+1);
    p2=blen(p+1+p1);
    r2=blen(r+1+r1);
    return ((p1==r1 && p2==r2 && equal(p+1,r+1) && equal(p+1+p1,r+1+r1)) ||
            (p1==r2 && p2==r1 && equal(p+1,r+1+r1) && equal(p+1+p1,r+1)));
  }
  if(p[0]==3) {
    p1=blen(p+1);
    r1=blen(r+1);
    p2=blen(p+1+p1);
    r2=blen(r+1+r1);
    p3=blen(p+1+p1+p2);
    r3=blen(r+1+r1+r2); 

    return ((p1==r1 && p2==r2 && p3==r3 && equal(p+1,r+1) && equal(p+1+p1,r+1+r1) && equal(p+1+p1+p2,r+1+r1+r2)) ||
            (p1==r1 && p2==r3 && p3==r2 && equal(p+1,r+1) && equal(p+1+p1,r+1+r1+r2) && equal(p+1+p1+p2,r+1+r1)) ||
            (p1==r2 && p2==r1 && p3==r3 && equal(p+1,r+1+r1) && equal(p+1+p1,r+1) && equal(p+1+p1+p2,r+1+r1+r2)) ||
            (p1==r2 && p2==r3 && p3==r1 && equal(p+1,r+1+r1) && equal(p+1+p1,r+1+r1+r2) && equal(p+1+p1+p2,r+1)) ||
            (p1==r3 && p2==r2 && p3==r1 && equal(p+1,r+1+r1+r2) && equal(p+1+p1,r+1+r1) && equal(p+1+p1+p2,r+1)) ||
            (p1==r3 && p2==r1 && p3==r2 && equal(p+1,r+1+r1+r2) && equal(p+1+p1,r+1) && equal(p+1+p1+p2,r+1+r1)));
  }
  else {printf("ERROR"); exit(0);}
}
 

int str[20][100000][60], nstr[20];

int main() {
  int N,n,b,k,i;
  printf("Antal kolatomer ?");
  scanf("%d", &N);
  nstr[0]=1;
  str[0][0][0]=0;
  for(n=1;n<N;n++) {
    nstr[n]=0;
    for(b=0;b<nstr[n-1];b++) {
      for(k=0;k<n;k++) if(str[n-1][b][k]<3) {
          for(i=0;i<k;i++) str[n][nstr[n]][i]=str[n-1][b][i];
          str[n][nstr[n]][k]=str[n-1][b][k]+1;
          str[n][nstr[n]][k+1]=0;
          for(i=k+1;i<n;i++) str[n][nstr[n]][i+1]=str[n-1][b][i];
          for(i=0;i<nstr[n];i++) if(equal(str[n][i],str[n][nstr[n]])) break;
          if(i==nstr[n]) {
            nstr[n]++;
          }
        }
    }
  }
   printf("Antal isomerer: %d\n", nstr[N-1]);
   return 0;
}

Här är en lösning i Java, som kanske inte är optimal (en del tunga operationer), men klarar ändå värsta testfallet (N=11) på 1 sek. (Man kan säkert göra en del optimeringar i koden också, men vem orkar göra det på en sådan här dryg uppgift.) Kruxet med den här uppgiften är helt klart att beskriva molekylerna. Krux nummer två är att tänka sig bort från en värd i 2D, det spelar ingen roll i vilken ordning "kolförgreningarna" sitter; om vi har en isomer som först förgrenar sig i tre grenar och sedan bara fortsätter på den vänstra grenen på ett unikt sätt, så är denna isomer den samma som en annan isomer där den unika fortsättningen sker på mittengrenen i stället. (Se isomererna i bilden för N=5. Betrakta isomeren längst ned till höger i rutan. Tänker man i 2D, så skulle vi helt klart få en ny isomer om vi bytte plats på grenen i mitten och den vänstra grenen, men i 3D skulle de båda isomererna anses lika. Det är lätt att tro att två isomerer bara är lika om de är symmetriska kring mitten, men så är inte fallet.)

När man väl har begripit vilka isomerer som ska anses som olika, så måste man beskriva isomererna på ett lämpligt sätt. I denna lösning görs det genom att beskriva hur många "barn" en viss kolatom har (hur många kolatomer den förgrenar sig till bortsett från kolatomen som förgrenar sig till den (den som ligger närmare OH-gruppen)), och sedan hur många barn barnen har i någon ordning. Med denna beskrivning klar måste man ta hänsyn till att det inte ska spela någon roll i vilken ordning man sedan väljer att beskriva barnen, vilket kan göras på 6 olika sätt; (V = vänster, M = mitten, H = höger) VMH, MHV, HVM, HMV, MVH, VHM. (Tänker man i 2D som jag först gjorde, så tror man att det bara är beskrivningarna VMH och HMV som behövs.) Dock kompliceras det ytterligare av att man inte nödvändigtvis behöver hålla sig strikt till en beskrivning hela tiden under en beskrivning av en isomer. Först kan man beskriva barnen i ordningen VMH, men Vs barn kan t.ex. beskrivas i ordningen HVM och Ms barn i ordningen VHM, och deras barn kan i sin tur beskrivas i godtycklig ordning.

Med beskrivningarna klara kan man bestämma att två isomerer ska anses lika när de har samma huvudbeskrivning. Det fiffiga är att i och med att två isomerer som är lika kommer att ha samma mängd beskrivningar (har de någon beskrivning lika kommer de att ha alla lika), så kan vi helt enkelt konsekvent bara välja den minsta beskrivningen hela tiden, för om de är lika så ska de också ha samma minsta beskrivning. På så sätt behöver vi bara hålla reda på en beskrivning per isomer, vilket är skonsamt för minnet. Och det än mer fiffiga är att detta inte bara gäller för isomeren i stort, utan också för delgrenarna i isomeren, så vi behöver egentligen inte generera alla beskrivningar, utan bara de minsta, så om vi har två delgrenar så behöver vi bara skapa kombinationer av deras respektive minsta beskrivningar, vilket snabbar upp programmet avsevärt. Sedan lagrar vi isomererna i en hashtabell så att dubbletter automatiskt försvinner. Först tar vi fram alla isomerer med 1 kolatom, sedan 2, sedan 3, osv och utifrån varje isomer skapar vi alla möjliga nya isomerer genom att lägga till en kolatom (tills vi fått så många vi önskar).

(Observera att det är viktigt att vi i Isomer-klassen överlagrar metoderna hashCode() och equals(). Skulle vi inte göra det, så skulle vi få heltokiga svar, eftersom båda metoderna i sådant fall skulle utgå ifrån minnesreferensen. För att objekt av ens egna klasser med fördel ska kunna stoppas in i en hashtabell måste man implementera metoderna så att följande gäller x.equals(y) => x.hashCode()==y.hashCode(); för om två objekt inte har samma hashvärde, så kan de inte hitta varandra och då kan man inte upptäcka att de är lika. Dvs eliminationen av dubbletter skulle inte fungera.)

import java.util.*;

public class Isomerer
{
	public static void main(String [] klein)
	{
		Scanner scan = new Scanner(System.in);

		System.out.print("Antal kolatomer: ");
		int atomer = scan.nextInt();

		//Sparar alla våra olika isomerer.
		HashSet <Isomer> isomerer = new HashSet <Isomer> ();

		//Lägger till och skapar grundisomeren (metanol).
		isomerer.add(new Isomer());

		//Då har vi använt oss av en kolatom.
		atomer--;

		//Så länge vi har kolatomer kvar, fortsätter vi att generera isomerer.
		for( ; atomer>0; atomer--)
		{
			//Sparar de isomerer vi genererar för tillfället.
			HashSet <Isomer> tmp = new HashSet <Isomer> ();

			//Så länge vi har isomerer kvar att behandla.
			for(Isomer i : isomerer)
			{
				//Genererar alla isomerar som kan genereras från denna isomer
				//genom att blott lägga till en enstaka kolatom.
				tmp.addAll(i.add());
			}

			//Lägger till de nya.
			isomerer = tmp;
		}

		//Skriver ut svaret.
		System.out.println("Antal isomerer: " + isomerer.size());
	}

	//En klass som beskriver en Isomer. En Isomer består internt av Kolatomer.
	//Klassen har en trädliknande struktur där isomeren blott refererar till "roten"
	//av isomeren (kolatomen ansluten till OH-gruppen) och kolatomerna i sin tur referarar
	//vidare till andra kolatomer i en riktning. En kolatom känner alltså inte till vem
	//den har som förälder. Man kan säga att vardera kolatom har tre barn, där några
	//(även alla) kan vara null.
	private static class Isomer
	{
		//Sparar alla olika sätt man kan lägga till en ny kolatom på.
		private LinkedList <String> toAdd = new LinkedList <String> ();

		//Huvudbeskrivningen (den "minsta" beskrivningen) av isomeren.
		String description = "";

		//Kolatomen som är den "första", den som är ansluten till OH-gruppen.
		KolAtom root;

		//Skapar en ny isomer med en kolatom (metanol).
		public Isomer()
		{
			root = new KolAtom();
		}

		//Skapar en ny isomer med den givna kolatomen (som kan ha andra kolatomer
		//redan anslutna till sig).
		private Isomer(KolAtom root)
		{
			this.root = root;
		}

		//Lägger till en ny kolatom till denna isomer.
		//Retunerar alla olika isomerer som kan uppstå av att man lägger till
		//en kolatom.
		public HashSet<Isomer> add()
		{
			//De olika isomerena som kan uppstå.
			HashSet <Isomer> versions = new HashSet <Isomer> ();

			//Hittar alla sätt man kan lägga till en kolatom på.
			add(root, "");

			//Lägger till en kolatom på alla dessa sätt.
			while(!toAdd.isEmpty())
			{
				//Plockar ut ett sätt att lägga till på.
				String way = toAdd.removeFirst();

				//Skapar en kopia av denna isomer.
				Isomer copy = new Isomer(root.clone());

				//Lägger till en ny kolatom på det specifierade sättet
				//på kopian av isomeren.
				copy.add(copy.root, way, 0);

				//Uppdaterar/skapar alla beskrivningar för denna isomer.
				copy.update();

				//Lägger till den i gänget som ska retuneras.
				versions.add(copy);
			}

			return versions;
		}

		//Tar reda på hur man kan lägga till en kolatom.
		//1 betyder att gå åt vänster.
		//2 betyder att gå åt mitten.
		//3 betyder att gå åt höger.
		//root: Kolatomen vi för tillfället är på.
		//way: Vägen vi har gått för att komma dit.
		private void add(KolAtom root, String way)
		{
			//Om kolatomen är null, så har vi funnit en plats
			//där vi kan lägga till en ny kolatom.
			if(root==null)
			{
				//Sparar vägen dit.
				toAdd.add(way);

				return;
			}

			//Går åt alla håll.
			add(root.left, way+1);
			add(root.center, way+2);
			add(root.right, way+3);
		}

		//Lägger till en ny kolatom på platsen som ges av way.
		//root: Kolatomen vi för tillfället är på.
		//way: Beskrivningen av vägen till platsen där vi ska lägga till atomen.
		//i: Håller reda på vilket index i way som talar om vilken väg vi ska
		//ta vid denna kolatom.
		private KolAtom add(KolAtom root, String way, int i)
		{
			//Om kolatomen är null, så har vi hittat platsen där vi ska lägga till en ny.
			//Tricket är också att den nya atomen kommer att retuneras och tilldelas till
			//referensen till den i föräldern.
			if(root==null) return new KolAtom();

			//Går vänster, mitten eller höger.
			if(way.charAt(i)=='1')
			{
				root.left = add(root.left, way, i+1);
			}
			else if(way.charAt(i)=='2')
			{
				root.center = add(root.center, way, i+1);
			}
			else
			{
				root.right = add(root.right, way, i+1);
			}

			//Oftast så ska ju referensen förbli den samma, så den retuneras till sig själv.
			return root;
		}

		//Uppdaterar beskrivningen för denna isomer. T.ex. efter en kolatom lagts till.
		public void update()
		{
			//Genererar det "minsta" sätt isomeren kan beskrivas på.
			description = root.update(root);
		}

		//Två isomerer är lika om de har samma beskrivning.
		public boolean equals(Object ob)
		{
			return ob instanceof Isomer && description.equals(((Isomer)ob).description);
		}

		//Hashvärdet avgörs av beskrivningen. På så sätt garanteras två
		//isomerer som är lika att ha samma hashvärde.
		public int hashCode()
		{
			return description.hashCode();
		}

		/********************************/
		//Klass som beskriver en kolatom.
		private static class KolAtom
		{
			//Länken till tre andra kolatomer.
			KolAtom left, center, right;

			public KolAtom()
			{
				//Skapar en kolatom.
			}

			//Retunerar en absolut kopia av denna kolatom inkl. kolatomerna den refererar till.
			public KolAtom clone()
			{
				try
				{
					KolAtom klon = (KolAtom)super.clone();
					if (left != null) klon.left = left.clone();
					if (center != null) klon.center = center.clone();
					if (right != null) klon.right = right.clone();

					return klon;
				}
				catch (CloneNotSupportedException e)
				{
					return null;
				}
			}

			//Genererar den minsta beskrivningen för denna kolatom tillsammans med de kolatomer den
			//refererar till. Beskrivningen byggs upp av en sträng, där det först står hur många
			//barn en kolatom har och sedan hur många barn dess olika barn har osv. Problemet är
			//i vilken ordning man ska ta barnens beskrivningar. För att vi ska kunna avgöra om
			//två isomerer som är lika men är uppbyggda olika är lika, så måste vi generera alla
			//möjliga beskrivningar. Först står hur många barn en kolatom har, detta gäller för
			//alla beskrviningar, men sedan kan det förekomma versioner där högerbeskrivningen
			//kommer före mittenbeskrivningen och vice versa (6 olika möjligheter), observera också
			//att ordningen av beskrivningarna för en nivå i "trädet"/isomeren inte behöver gälla för
			//en annan nivå eller del av isomeren; detta gör genereringen av beskrivningar än mer komplicerad.
			//
			//Hur som helst retunerar vi den minsta beskrivningen (som ges av hur strängar jämförs).
			//Tricket i denna rekursiva metod är att den minsta beskrivningen som går att få består av de minsta
			//delbeskrivningarna som går att få. Och den minsta beskrivningen är givetvis den minsta kombinationen
			//av de minsta delbeskrivningarna.
			public String update(KolAtom root)
			{
				//Den icke existerande kolatomen beskrivs logiskt nog som ingenting.
				if(root==null) return "";

				//Hur många barn denna kolatom har.
				int barn = 0;

				//Räknar barn.
				if(root.left!=null) barn++;
				if(root.center!=null) barn++;
				if(root.right!=null) barn++;

				//Beskrivningen för en änd-kolatom är blott 0 då den har 0 barn.
				if(barn==0)
				{
					return "0";
				}

				//Hämtar minsta beskrivningen för respektive barn.
				String combs1 = update(root.left);
				String combs2 = update(root.center);
				String combs3 = update(root.right);

				//Genererar den minsta beskrivningen som går att kombinera ihop.
				String comb = comb(""+barn, combs1, combs2, combs3);

				//Retunerar den minsta beskrivningen.
				return comb;
			}

			//Retunerar den minsta kombinationen/permutationen som går att få
			//där strängen barn ska stå först.
			private String comb(String barn, String c1, String c2, String c3)
			{
				//Det finns 6 olika sätt att kombinera ihop dem på.
				//Detta är ett sätt.
				String comb = barn+c1+c2+c3;

				//Detta är ett annat sätt.
				String tmp = barn+c2+c3+c1;

				//Om det andra sättet är mindre än det första, så sparar
				//vi det andra sättet och gör oss av med det första.
				//Annars tvärtom.
				comb = (tmp.compareTo(comb)<0) ? tmp : comb;

				//Upprepar samma procedur för de 4 övriga sätten.
				tmp = barn+c3+c1+c2;
				comb = (tmp.compareTo(comb)<0) ? tmp : comb;
				tmp = barn+c3+c2+c1;
				comb = (tmp.compareTo(comb)<0) ? tmp : comb;
				tmp = barn+c2+c1+c3;
				comb = (tmp.compareTo(comb)<0) ? tmp : comb;
				tmp = barn+c1+c3+c2;
				comb = (tmp.compareTo(comb)<0) ? tmp : comb;

				//Retunerar den minsta beskrivningen.
				return comb;
			}
		}
	}
}

Referenser[redigera]

↑ Cayley, A. (1874) On the Mathematical Theory of Isomers. Philos. Mag., 47, 444
↑ Paton, R. S., Goodman, J. M. (2007) Exploration of the Accessible Chemical Space of Acyclic Alkanes. J. Chem. Inf. Model., 47, 2124-2132

[first-1] Cayley, A. (1874) On the Mathematical Theory of Isomers. Philos. Mag., 47, 444

[rek-2] Paton, R. S., Goodman, J. M. (2007) Exploration of the Accessible Chemical Space of Acyclic Alkanes. J. Chem. Inf. Model., 47, 2124-2132

[1]

[2]