Výpočet determinantu v pythnu

28.2.2008 21:19 | Přečteno: 4744× | Další věci

Dělám do školy program na výpočet libovolného determinantu. Když jsem po týdnu vymyslel jak na to, tak mi to někde blbne. Nacpu tomu list (př [[1, 2, 3][4, 8, 7][8, 9, 6]]) vnořený list je řádek)

def detr(det, i): #Tohle upravuje determinant dělá z něj sub determinant
    K = det
    del K[0]
    for y in xrange(len(K)):
        del K[y][i]
    return K

def determinant(det, stupen):
    if stupen == 2:
        d = det[0][0]*det[1][1] - det[0][1]*det[1][0]
    else: 
        d = 0
        for i in xrange(stupen):
            M = detr(det, i)
            print stupen
            d = d + (-1)**(i-1)*det[0][i]*determinant(M, (stupen -1) ) 
    return d

Až po první determinant 2. stupně to jede dobře, jenže když to má udělat další determinant 2. stupně, tak se to K zase zkrátí a já už prostě nevím, jak udělat to zkracovaní. Předem moc dík za každou radu. Vím, že existuje numpy, ale v zadání je zpracujte algoritmus. Kdyby měl někdo jiný nápad na řeší determinantů přivítal bych to.        

Komentáře

[vencour]

Momentálně si vzpomenu jen na výpočet rozkladem a to by se snad dělalo rekurzí, když neznáš předem hodnost matice?

[Saljack]

Vždyť to je rekurze. A ten rozklad to je snad taky je to to, že se vyškrtne jedna řádka a sloupec vznikne ti determinant ten vynásobíš (-1)**(i-1)*a(11) pak vyškrtneš změníš přeškrtnutý sloupec a zase uděláš to samí co před tím

[vencour]

Ok, do pythonu jsem se ještě nedostal. Ten součet permutací, sudejch a lichejch sloupců a řádků asi zohledňuješ?

[Saljack]

Permutaci neumím, tak daleko nejsem

[vencour]

Mluvim o tom, že na řádku 1, skoupci 1, je determinant kladnej, na 1/2 je zápornejm na 2/2 kladnej ...

[Saljack]

jj to je to (-1)**(i-1)*det[0][i] ty ** je mocnina

[Ignor]

Uá. Ty jsi teda hustilka. Aneb máme rádi faktoriál.

[Ignor]

Tak ať nejsem karpatnej jako ten hrad v Tejemnech: http://en.wikipedia.org/wiki/Determinant#Algorithmic_implementation

[Saljack]

Já ten výpočet mám dobře (alespoň doufám), akorát nevim jak do toho zadělat to vyškrtávání sloupců a řádků

[Jirka P]

Python neznám, ale podle mně jde o to, že když napíšeš

y = [1, 2, 3]
x = y
del x[0]
if y[0] == 2 :
   print "true"

tak to vypíše true, protože přiřazení je přiřazení referencí jako v Javě. Dá se to spravit něčím jako x=y[:]. Doporučuju přečíst si nějakou učebnici Pythonu, tohle by měl být základ.

[Saljack]

Díky díky moc díky měl jsi pravdu já o tom vím ale úplně jsem na to zapoměl

[petr]

zapoměl

Tak k té příručce Pythonu ještě naštuduj mluvnici, abys na to příště nezapomněl

[Martin Tůma]

Možná že to máš "dobře" (i když když to nefunguje tak asi těžko...), ale prakticky je to zcela nepoužitelný. Protože u algoritmu se složitostí O(n!) končíš na n v řádu jednotek.

Takže ten kód rovnou zahoď a napiš si prachobyčejnou Gaussovu eliminaci. Vynásobenim prvků na diagonále pak získáš determinant. Složitost O(n^3).

[Josef Kufner]

Napiš si mapovací funkci. Té funkci dáš parametry řádek (r), sloupec (s), vyškrtnutý řádek (vr) a vyškrtnutý sloupec (vs). A ono ti to vrátí prvek v matici o požadovaných souřadnicích, kde přepočítané souřadnice [pr, ps] budou vypadat nějak takto (psáno v C):

pr = r < vr ? r : r + 1;
ps = s < vs ? s : s + 1;

Ten zbytek už je jednoduchý.

[Ignor]

A v dalším kroku (dva vyškrtnuté sloupce)?

Do matice přistupuj nepřímo přes vektor indexů. Takový vektor ti stačí jeden - pro sloupce nebo řádky, ale musíš rozvíjet vždy podle prvního řádku nebo prvního sloupce (v kroku k, k-ého sloupce _původní_ matice).

V indexovém vektoru můžeš buď vyšktrtávat (měnit ho, měnit jeho délku), nebo jen prohazovat prvky v(1)..v(k) jsou nevyškrtané sloupce, v(k+1)..v(n) vyškrtané, musíš si teda pamatovat k. Vyškrtnutí sloupce j: a = v(j); v(j) = v(k); v(k) = a; k = k - 1;

Ale fakt radši tu gaussovku...

[Saljack]

Já nechápu co tady všichni máte s gaussovkou, to by bylo úplně na nic tím nevypočtu ten determinant. Gaussovku budu psát pro matice.

[Ignor]

Kdybyses podíval na ten odkaz, co jsem ti dal...

BTW až nakonec budeš tu G. elim. psát, tak nezapomínej, co s determinantem udělá vynásobení řádku koeficientem atd.

A když něco nechápeš (jak sám říkáš), tak neříkej, že je to na nic, na tohle jsem háklivej

[Saljack]

ale vždyť to je ono jestli jsi myslel tohle

[Ignor]

Zkus na minutu předpokládat, že nejsem úplnej debil. Teď se pokus znovu navštívit onen odkaz Ujisti se, že koukáš na část 9 Algorithmic implementation hesla Determinant na americké Wikipedii. Pak se pozvolna vrať k tomu, že jsem debil

[Saljack]

Já neříkám, že jsi debil. Spíš já jsem debil. Už jsem to pochopil. Ale měl jsem taky pravdu gaussovou eliminační metodou jenom tu matici upravíš na trojúhelníkoví tvar a z toho potom vycházíš že det(M) = det(LU) = det(L)det(U). Jenže zrovna tohle díky zadání nejde použít . Moc dík že tady se mnou ztrácíš čas

[Martin Tůma]

Už jsem to pochopil.

Nepochopil. Zkus to tady, třeba tomu budeš rozumět v češtině líp.

Pokuď v GEM používáš pouze úpravu typu přičtení násobku řádku k jinému řádku jak je obvyklé, získáš determinant původní matice pouhým vynásobením prvků na diagonále trojúhelníkové matice vzniklé z matice původní pomocí GEM.

Nevím jaký je přesně zadání, ale pokuď v něm neni, že se MUSÍ použít výpočet pomocí subdeterminantů, tak výše uvedenj postup použít rozhodně lze. Naopak, pokuď v zadání není dáno, že musíš použít výpočet pomocí subdeterminantů, měl by tě s tím co tady "patláš" zadavatel okamžitě vyhodit. Protože s někým, kdo použije algoritmus s O(n!) tam, kde existuje jednodušší algoritmus s O(n^3) se nic jinýho udělat nedá...

[tsLnox]

1!

Nemám rád faktoriál.
Nemám faktoriál rád.
Rád nemám faktoriál.
Rád faktoriál nemám.
Faktoriál nemám rád.
Faktoriál rád nemám.

xD

(nemám z vlastní hlavy...)

[freshmouse]

Najdi si Petra Olšáka a jeho skriptum Úvod do algebry, zejména lineární. Je ke stažení na webu (v tištěné verzi tam jsou i polynomy, ale možná jsou už i na webu).

... nebudu tě trápit, je to tady.

[freshmouse]

Olšák tam má i návod na řešení soustav rovnic strojově (str. 60), možná by ti to mohlo pomoct. A nezapomeň si přečíst, jak se mění determinant po použití Gaus. elim. m. (str. 51). Možná to nepochopíš z té definice, ale pak je tam i příklad (př. 4.23 + pozn. 4.22), a z toho by ti to mělo být jasné.

[KS]

Tyjo to jsem dělal ještě na střední ve druháku. Můžu dát jedině kód v C#, ale nijak za něj neručím, v té době jsem ještě ani nevěděl, co je to lineární algebra a C# jsem taktéž moc neovládal. Takže nepište že jsem vůl a tak.

using System;

namespace Calc
{

	public class Matrix
	{
		protected double[,] MatrixData;
		protected int my, nx;
		
		//
		//INDEXER
		//
		//for access to mambers
		//arrays are indexing from 0! but matrix from 1
		public double this [int m, int n]
		{
			get
			{
				return MatrixData[m-1, n-1];
			}
			set
			{
				MatrixData[m-1, n-1] = value;
			}
		}
	
		//
		//CONSTRUCTORS
		//
		//constructor with size of MatrixData
		public Matrix(int m, int n)
		{
			//m+1 and n+1 for counting from 1 (matrix is a11 .. amn)
			MatrixData = new double[m, n];
			my = m;
			nx = n;
		}

		//matrix from array and size
		public Matrix(double[,] MatrixArray, int m, int n)
		{
			MatrixData = new double[m, n];

			for(int x = 0; x <= m-1; x++)
				for(int y = 0; y <= n-1; y++)
					MatrixData[x, y] = MatrixArray[x, y];

			my = m;
			nx = n;
		}

		//new matrix from existing matrix
		public Matrix(Matrix ExistingMatrix)
		{
			MatrixData = ExistingMatrix.GetMatrix();
			my = ExistingMatrix.GetM();
			nx = ExistingMatrix.GetN();
		}

		//
		//OPERATORS
		//
		//addition Matrix+Matrix
		public static Matrix operator+(Matrix A, Matrix B)
		{
			//size of matrix must be the same as A and B
			Matrix C = new Matrix(A);

			for(int m = 1; m <= A.GetM(); m++)
				for(int n = 1; n <= A.GetN(); n++)
					C[m, n] = A[m, n] + B[m, n];

			return C;
		}

		//multiply Matrix*Matrix
		public static Matrix operator*(Matrix A, Matrix B)
		{
			Matrix C = new Matrix(A.GetM(), B.GetN());
			double sum;

			for(int m = 1; m <= A.GetM(); m++)
				for(int n = 1; n <= B.GetN(); n++)
				{
					sum = 0;

					for(int i = 1; i <= A.GetN(); i++)
						sum += A[m, i] * B[i, n];

					C[m, n] = sum;
				}

			return C;
		}

		//multiply double*Matrix
		public static Matrix operator*(double a, Matrix B)
		{
			Matrix C = new Matrix(B);

			for(int m = 1; m <= B.GetM(); m++)
				for(int n = 1; n <= B.GetN(); n++)
					C[m, n] = a * B[m, n];

			return C;
		}

		//multiply Matrix*double
		public static Matrix operator*(Matrix A, double b)
		{
			return b * A;
		}

		//compare
		public static bool operator==(Matrix A, Matrix B)
		{
			return A.Equals(B);
		}

		//compare
		public static bool operator!=(Matrix A, Matrix B)
		{
			return !A.Equals(B);
		}

		//invert
		public static Matrix operator~(Matrix X)
		{
			return X.Invert();
		}
		
		//
		//GETS
		//
		//returns array with matrix
		public double[,] GetMatrix()
		{
			return (double[,])MatrixData.Clone();
		}

		//returns m - y size
		public int GetM()
		{
			return my;
		}

		//returns n - x size
		public int GetN()
		{
			return nx;
		}

		public override bool Equals(object obj)
		{
			if(obj is Matrix)
			{
				Matrix A = (Matrix)obj;

				for(int m = 1; m <= this.GetM(); m++)
					for(int n = 1; n <= this.GetN(); n++)
						if(this[m, n] == A[m, n])
							return true;

				return false;
			}
			else
				return false;
		}

		//
		//COUNTING
		//
		//second determinant
		public double Determinant2()
		{
			return this[1, 1] * this[2, 2] -
				this[1, 2] * this[2, 1];
		}

		//third determinant
		public double Determinant3()
		{
			return  this[1, 1] * this[2, 2] * this[3, 3] -
				this[1, 1] * this[2, 3] * this[3, 2] +
				this[1, 2] * this[2, 3] * this[3, 1] -
				this[1, 2] * this[2, 1] * this[3, 3] +
				this[1, 3] * this[2, 1] * this[3, 2] -
				this[1, 3] * this[2, 2] * this[3, 1];
		}

		//n-scale determinant
		public double Determinant(int n)
		{
			switch(n)
			{
				case 2: 
					return Determinant2();
				case 3:
					return Determinant3();
				default:
					double result = 0;
					for(int i = 1; i <= this.GetN(); i++)
						result += this[1, i] * this.Addition(1, i);
					return result;
			}
		}

		//addition of mamber
		public double Addition(int m, int n)
		{
			Matrix X = new Matrix(this.GetM()-1, this.GetM()-1);
			int i = 1, j;

			for(int x = 1; x <= this.GetM(); x++)
			{
				j = 1;
				for(int y = 1; y <= this.GetM(); y++)
				{
					if(m != x && n != y)
					{
						X[i, j] = this[x, y];
						j++;
					}
				}
				if(x != m)
					i++;
			}

			return System.Math.Pow(-1, m+n) * X.Determinant(this.GetM()-1);
		}

		//invert matrix
		public Matrix Invert()
		{
			Matrix X = new Matrix(this.GetM(), this.GetM()); 

			for(int m = 1; m <= this.GetM(); m++)
				for(int n = 1; n <= this.GetN(); n++)
					X[n, m] = this.Addition(m, n);

			return (1/this.Determinant(this.GetM())) * X;
		}

		//rotate matrix right
		public void RotateRight()
		{
			Matrix X = new Matrix(this.GetN(), this.GetM());

			for(int m = 1; m <= this.GetM(); m++)
				for(int n = 1; n <= this.GetN(); n++)
					X[n, this.GetM() - m + 1] = this[m, n];

			MatrixData = X.GetMatrix();
			my = X.GetM();
			nx = X.GetN();
		}

		public Matrix Transpone()
		{
			Matrix X = new Matrix(this.GetN(), this.GetM());

			for(int m = 1; m <= this.GetM(); m++)
				for(int n = 1; n <= this.GetN(); n++)
					X[n, m] = this[m, n];

			return X;
		}

		public override string ToString()
		{
			string s = "";

			for(int m = 1; m <= this.GetM(); m++)
			{
				s += "(";
				for(int n = 1; n < this.GetN(); n++)
					s += this[m, n].ToString() + "; ";

				s += this[m, this.GetN()].ToString() + ")\n";
			}

			return s;
		}

		public object Clone()
		{
			return new Matrix(this);
		}
	}
}

[petr]

No, pořád mi ten zdroják v C# přijde čitelnější než ten v Pythonu. A to mám cé žbrdlinku rád asi tak jako pana Sponku! Nevím, ale mít v hlavičce funkce jen identifikátor bez udaného typu -- a ještě bez komentářů ??? To na mě opravdu moc není

[KS]

Syntaxe C# mi přijde hodně intuitivní. Když zapomenu nějakou syntaxi, tak než abych to hledal v dokumentaci, tak tam prostě zkusím něco prásknout, jak by to asi mohlo být, a většinou se to přeloží. Teď jsem programoval v C++/CLI a vůbec to tak intuitivně nefungovalo. V jednom kuse jsem študoval dokumentaci. Tahle intuitivnost mi přijde důležitá.

[freshmouse]

Čitelnost C# mi přijde stejná jako čitelnost Javy, Pythonu a dalších jazyků. Rozdíl oproti Javě je třeba ten, že "import" nahradili "using", "package" nahradili "namespace" atd. Třeba cykly, podmínky a další věci jsou ve většině (konvenčních ) jazyků stejné, a jestli použiješ ".", nebo "->", to je jen o zvyku. Stejně tak ukončení příkazu: jestli středníkem, koncem řádky -- to už je jedno.

Nejsem nějaký velký programátor a přestože znám jazyků jen pár, jsem schopen číst většinu ostatních (počítají se jen ty "normální", většinou se syntaxí odvozenou z C, takže ne Perl a pak různé LISPy atd. nic ve zlém...).

[KS]

Jasně, já jsem měl spíš na mysli intuitivnost syntaxe. Když si přesně nepamatuješ syntaxi, ale dokážeš ji nějak odhadnout z toho co víš aniž bys musel hledat někde v knihách. Třeba u toho C++/CLI je hodně klíčových slov, třeba generic a typename u generických tříd a to si prostě musíš naštudovat.

[freshmouse]

Spíš myslíš intuitivnost názvů standardních balíků, tříd, metod atd. než intuitivnost syntaxe.

[KS]

No tak jo. Když víš, co si myslím .

[Kyosuke]

Tak v tomhle směru jsem zjistil, že myšlenky mi čte autor Gauche Scheme. Já myslím, že každý má tu intuici seřízenou trošku jinak.

[Kvakor]

Pokud vim, tak citelnost byla dokonce z jednim z cilu pri navrhu C#. Kdezto Python je jen "hezci Perl" (doufam, ze jsem zadne Pythonisty neurazil) a Perl rozhodne nepatri mezi jazyky, ktere jsou snadno citelne (dokonce existoval i Obfuscated Perl Contest, coz byla obdoba The International Obfuscated C Code Contest).

BTW: žbrdlinka je nadherne slovo

[freshmouse]

Ach jo.

A co bys jako chtěl, aby ti tvůrci C# (Javy, Pythonu...) řekli? Že cílem návrhu byla nečitelnost?

Neznám Python ani Perl, ale Python opravdu není "jen hezčí Perl".

Koneckonců, on ten Perl má svoje kouzlo, pokud na jeho hru přistoupíš. Pak se vyjadřuješ úsporně a rychle. Jeho jediným problémem je, že se v něm asi dá těžko dělat v týmu...

V Perlu můžeš psát i "běžně", ale to pak zrácíš jeho výhody při vyjadřování.

[MJ]

Jeho jediným problémem je, že se v něm asi dá těžko dělat v týmu...

Řekl bych, že existence CPANu to dosti přesvědčivě vyvrací

Neznám Python ani Perl, ale Python opravdu není "jen hezčí Perl".

To zrovna nezní jako kvalifikovaný názor

[freshmouse]

Mno, CPAN není nic v týmu, tam může přispět každý -- svým oblíbeným stylem zápisu programu... Že neznám Python a Perl znamenalo, že je neznám dokonale, neprogramuju v nich, ale něco málo bych o nich věděl.

[happy barney]

hmm, to znie, ako keby bol perl nejaký mystický jazyk

to nie, on je len obeťou masívnej ohováracej kampane tých, čo mu nie sú schopní porozumieť

[thingie]

Tak on by se ten pythonový dal napsat hezčejš. Ta C# je IMHO jenom vleklá onanie s for smyčkou :o)

[Petr "Glubo" Sýkora]

Jsi vůl ne kvuli obsahu toho zdrojáku, ale kvůli tomu, že jsi ho nepřidal jako přílohu ke příspěvku. A prosím nebrat vážně tvrzení, že jsi vůl, nebýt toho, žes to sám použil, nenazval bych tě tak (proto ten smajlík)

[MJ]

Pěkné, ale než to doběhne na matici 30x30, tak bych se s tím raději nechlubil

Používá to totiž úplně stejný naprosto nešikovný algoritmus jako ten, o kterém je řeč.

[Kyosuke]

Hmm, osobně bych to popsal slovy "pro samá public static nevidím algoritmus".

[KOLiK]

No pokud to máš jako cvičení, tak budiž :) Ale pyton už pěkné knihovny pro matematiku má, a proto není potřeba už takové elementární věci vymýšlet. Stačí použít numpy a scipy (která numpy obsahuje):

from numpy import matrix
from scipy.linalg import det

A = matrix([[1,2],[3,4]])
print det(A)

[KOLiK]

Aha, mměl jsem si přečíst příspěvěk pořádněji ;)

[CET]

Predem priznam, ze Python zatim jenom okukuju. Ale zkusil jsem si tu tvoji cast hodit do textaku a drobne doplnit a upravit a zkusit. Udelal jsem tohle

#!/usr/bin/python

def detr(det, i): 
  K = det
  del K[0]
  for y in xrange(len(K)):
    del K[y][i]
  return K

def determinant(det, stupen):
  print "Stupen1: " , stupen
  if stupen == 2:
    print "Stupen2a: " , stupen
    d = det[0][0]*det[1][1] - det[0][1]*det[1][0]
  else: 
    print "Stupen2b: " , stupen
    d = 0
    print "len det 0 " , len(det[0])
    for i in xrange(stupen):
      M = detr(det, i)
      print "det[0][",i,"]=",det [0][i]
      print "M=",M
      print "det=",det
      d = d + (-1)**(i-1)*det[0][i]*determinant(M, (stupen -1) ) 
  print "Subdeterminant pro " , stupen , " je " , d
  return d

mat=[[1, 2, 3],[4, 8, 7],[8, 9, 6]]
print determinant(mat,len(mat));

Bohuzel mi to ale hlasi chybu - tohle je cely vystup z toho programu.

Stupen1:  3
Stupen2b:  3
len det 0  3
det[0][ 0 ]= 8
M= [[8, 7], [9, 6]]
det= [[8, 7], [9, 6]]
Stupen1:  2
Stupen2a:  2
Subdeterminant pro  2  je  -15
det[0][ 1 ]=
Traceback (most recent call last):
  File "./det.py", line 30, in ?
    print determinant(mat,len(mat));
  File "./det.py", line 21, in determinant
    print "det[0][",i,"]=",det [0][i]
IndexError: list index out of range

Pokud jsem python dobre pochopil a podle vypisu interpretoval, tak python predava seznamy jako odkazy a nedela kopie. Cili, jakakoliv destruktivni zmena seznamu K (resp. parametru "det") ve funkci detr se projevi zmenou puvodniho seznamu "det" ve funkci determinant. Jak jsem psal, nejsem pythonista, takze ti neporadim, jak udelat kopii seznamu a predat do funkce k destrukci tak, aby puvodni seznam zustal, ale nejakej pythonista by to vedet mohl. Zkusim jeste mrknou na nejaky ucebnice, neco by se najit mohlo.

[CET]

Takze, protoze nejsem pythonista, ale perlic, tak jsem zaperlil a nechal jsem se nachytat i perlu. Samozrejme, ze pokud chci v perlu pouzit vicerozmerny pole, tak musim pouzit reference. Jinak bych musel mit jednorozmerny pole a to by se pak do funkce predavalo jako hodnota (zkopirovaly by se vsechny prvky). Kdyz mam ale dvojrozmerne pole, tak vnitrni pole museji byt jako anonymni array, cili reference. A jakmile zacnu mazat v referencich, tak je to na nic. (pokud je tu nejakej dalsi perlista a nevic profik, tak me dyztak oprav, jestli se daji vytvorit dvojrozmerny pole referencovany pomoci $a[x][y], aniz by prvny pole @a byly reference na pod-pole).

Takze jsem splodil nasledujici:

#!/usr/bin/perl

sub detr {
  my ($i,$k)=@_;
  my $r=[];
  for (my $y=1;$y<scalar(@{$k});$y++) {
    push @{$r},[@{$k->[$y]}];
    splice @{$r->[$y-1]},$i,1;
  }
  return $r;
}

sub determinant {
  my ($det)=@_;
  my $stupen=scalar(@{$det});
  if ($stupen==1) {
    return $det->[0][0];
  } elsif ($stupen==2) {
    return $det->[0][0]*$det->[1][1]-$det->[0][1]*$det->[1][0];
  } else {
    my $d=0;
    for (my $i=0;$i<$stupen;$i++) {
      my $m=detr($i,$det);
      $d+=(-1)**($i)*$det->[0][$i]*determinant($m);
    }
    return $d;
  }
}

$mat=[
  [1,2,3],
  [4,8,7],
  [8,9,6]
];

print determinant($mat) . "\n";

A mam k tomu par komentaru:

• Bacha na ty reference, jak jsem psal v predchozim postu, python predava pole (ted nevim, jestli jenom vnitrni pole, nebo i to hlavni) jako referenci, tak upravy se provadeji v tom hlevnim objektu, ktery tim padem nicis a proto mi to psalo chybu "list index out of range".
• dalsi vec je v tom vzorecku - mas tam chybu, mas tam (-1)**(i-1), ale spravne ma byt (-1)**i, protoze ten vzorecek je spravne (-1)**(radek+sloupec) takze zacinas 0+0, coz je nula a -1 na nultou je jedna kladna, takze je to dobre, ale tobe vychazi pro prvni hodnotu, kdy i==0, -1 na -1, coz je -1 a to je spatne.
• No a posledni je uz jenom rypnuti Pouzitim konkretniho jazyka se z tebe lepsi programator nestane. On si kazdej mysli, ze kdyz se nauci python, tak bude umet programovat. Neni to tak, je uplne jedno, jaky jazyk pouzijes, pokud ten jazyk neznas dobre, nechapes vnitrni funkci a neumis vytvorit spravne algoritmus, tak to za tebe zadny jazyk neudela. Mozna navic, kdybys tohle delal v Cecko (nemyslim C+ ani C#, ale klasicky Cecko), pripadne rovnou assembler (teda hlavne tim myslim tu praci s polema), tak bys chybu nasel zrejme hned. Bohuzel ten python oslepil (nebo mozna lidi, kteri ho presprilis vychaluji) a nevidis, co jsi skutecne naprogramoval. Takze klidne v pythonu pokracuj, ale nauc se ho poradne.

BTW: Vysledek u tyhle matice je -35, overeno v C-style arbitrary precision calculator (version 2.12.1.5).

[CET]

Tak jsem jeste koukal, jak by se v pythonu dalo zkopirovat pole, na ktery mam odkaz. Neco jako v perlu:

$a=[ [1,2,3] ];
$b=[ [@{$a->[0]}] ];

ono to sice vypada neprehledne, ale pokud clovek zna aspon zaklady perlu, tak to uvidi hned. $a je reference na anonymni array, ktery obsahuje jeden prvek a tim je odkaz na anonymni array se tremi prvky 1, 2 a 3. Na druhem radku pak vytvarim do promenne $b referenci na jiny anonymni array, do ktereho vlozim jeden prvek, ktery je odkaz na dalsi anonymni pole, do ktereho vlozim hodnoty - a ty se vlozit zpusobem takovym, ze se odkazu na 1.prvek (index 0) anonymniho pole odkazovaneho promennou $a $a->[0], tento prvek obsahuje odkaz na anonymni pole s prvky 1, 2, 3, takze to necham vyhodnotit jako obycejne pole @{$a->[0]}, coz bude prave to pole tech hodnot (uz to tedy nebude odkaz na to anonymni pole, ale konkretni -zkopirovane- hodnoty) a ty se mi vlozi do toho noveho anonymniho pole [ @{$a->[0]} ], ktere je ulozene jako prvni prvek celeho anonymniho pole odkazovaneho prom. $b = [ [ @{$a->[0]} ] ]

Kdyz se pokusim v pythonu udelat tohle

>>> a=[ [ 1,2,3] ]
>>> b=[]     
>>> c=[]
>>> b+=[a[0]]
>>> c.append(a[0])
>>> a
[[1, 2, 3]]
>>> b
[[1, 2, 3]]
>>> c
[[1, 2, 3]]
>>> del b[0][1]
>>> del c[0][1]
>>> a
[[1]]
>>> b
[[1]]
>>> c
[[1]]

vidite sami, ze oba zpusoby b+=... i c.append nevlozi konkretni hodnoty, ale referenci. V docu k pythonu to pisou, ze se kopiruje reference, ale nenasel jsem nic o tom, jak seznam vratit jako novy seznam hodnot.

[CET]

Hehe a jeste jedna vec. Ted me napadlo, ze tu rekurzi lze hnat az do extermu a netestovat ten stupen 2, ale nechat to az na stupen 1, kde se vraci pouze jedine cislo.

#!/usr/bin/perl

use Data::Dumper;

sub detr {
  my ($i,$k)=@_;
  my $r=[];
  for (my $y=1;$y<scalar(@{$k});$y++) {
    push @{$r},[@{$k->[$y]}];
    splice @{$r->[$y-1]},$i,1;
  }
  return $r;
}

sub determinant {
  my ($det)=@_;
  my $stupen=scalar(@{$det});
  if ($stupen==1) {
    return $det->[0][0];
  } else {
    my $d=0;
    for (my $i=0;$i<$stupen;$i++) {
      my $m=detr($i,$det);
      $d+=(-1)**($i)*$det->[0][$i]*determinant($m);
    }
    return $d;
  }
}

$mat=[
  [1,2,3],
  [4,8,7],
  [8,9,6]
];

print determinant($mat) . "\n";

[podlesh]

Nejsem vůbec pythonista (no, asi před čtyřmi roky jsem napsal takového démonka), ale nemohu se zbavit dojmu že hledáš problémy tam kde nejsou.

a=[[1], [2], [3], [4], [5, "a"]]
for i in range(0,len(a)):
     b=a[:i]+a[i+1:]
     print b

Jinak vytvoření kopie listu není problém:

b=list(a)

[mimi.vx]

kopie prommenne v pythonu: a=[1,2,3] b=a[:]

a nebo modul copy a funkce deepcopy:)

[CET]

jo, ten operator [i:j] slice jsem samozrejme neznal, ale vypada to, ze by to mohlo byt reseni. Takze jsem to trosku predelal a vypada to nasledovne:

varianta pro konecny stupen 2:

#!/usr/bin/python
# vim:set et:

def detr(det, i): 
  K = []
  for y in range(1,len(det)):
    K.append(det[y][:i]+det[y][i+1:])
  return K

def determinant(det, stupen):
  if stupen == 2:
    d = det[0][0]*det[1][1] - det[0][1]*det[1][0]
  else: 
    d = 0
    for i in xrange(stupen):
      M = detr(det, i)
      d = d + (-1)**i*det[0][i]*determinant(M, (stupen -1) ) 
  return d

mat=[[1, 2, 3],[4, 8, 7],[8, 9, 6]]
print determinant(mat,len(mat));

a pak jeste pokud bys chtel koncit az ve stupni 1

#!/usr/bin/python
# vim:set et:

def detr(det, i): 
  K = []
  for y in range(1,len(det)):
    K.append(det[y][:i]+det[y][i+1:])
  return K

def determinant(det, stupen):
  if stupen == 1:
    d = det[0][0]
  else: 
    d = 0
    for i in xrange(stupen):
      M = detr(det, i)
      d = d + (-1)**i*det[0][i]*determinant(M, (stupen -1) ) 
  return d

mat=[[1, 2, 3],[4, 8, 7],[8, 9, 6]]
print determinant(mat,len(mat));

No, vypada to, ze to chodi. Tak snad to pomuze.

[Jary]

Možná by nebylo od věci použít array, což je něco jako Cčkové jednotypové pole. Obalilo by se to nějakejma setter-getter metodama a kopie by pak byla jednoduchá a mělo by to bejt rychlé a nenažrané.

[Martin Tůma]

Příloha:

• determinant.png (3053 bytů)

Vždycky mě dostane do kolen, když script-kiddies hodlají "optimalizovat" svůj algoritmus s exponenciální složitostí pro problém, kde existuje řešení s polynomiální složitostí pomocí podobnejch metod...

Pro ty co by snad na tuhle diskuzi narazili někdy v budoucnu a chtěli odsud kopírovat některý z uvedených kódů v příloze uvádim obrázek, po jehož shlédnutí by si to měl rozmyslet i negramotný chovanec Guatemalského ústavu pro choromyslné...

[CET]

Ale jo, samozrejme, ze jsou lepsi zpusoby. Ja osobne bych to resil instalaci calcu a zadanim

mat A[X,X]=[ hodnoty matice]
det(A)

a mam vypocet rychle a presne (prepokladam, ze v takovym programu budou pouzity jedny z nejrychlejsich algoritmu).

Nicmene, hledani chyb v cizich programech je dobra skola programovani. Clovek u toho musi premyslet, co tim chtel autor vlastne rict a kde je teda chyba.

[freshmouse]

V čem jsi dělal ten graf?

[Käyttäjä 11133]

Gnuplot... IMHO

[Martin Tůma]

V gnuplotu.

set terminal png size 500,300
set output "determinant.png"
set xrange [0:8]
set grid
set xlabel "Velikost matice"
set ylabel "Čas výpočtu"
plot x**3 title "Gaussova eliminace", \
     sqrt(2*pi*x)*x**x * exp(-x) title "Metoda kofaktorů"

[Martin Tůma]

Příloha:

• det_logscale.png (3959 bytů)

Možná ještě jeden "údernější":

set terminal png size 500,300
set output "determinant.png"
set xrange [2:100]
set logscale y
set grid
set title "Porovnání výpočetní složitosti determinantu matice"
set xlabel "Velikost matice"
set ylabel "Čas výpočtu"
plot x**3 title "Gaussova eliminace", \
     sqrt(2*pi*x)*x**x * exp(-x) title "Metoda kofaktorů"

[Käyttäjä 11133]

Příloha:

• jeste-udernejsi.png (11442 bytů)

A ještě údernější

[Jary]

Tak to vypadá, že jsem si naběhl na smeč (Hmmm, nejsem si jist, že tuhle frázi chápu dobře) tím, že mě zaujal jen ten problém s deepcopy. Řešení Gaussovou eliminací, kterou jsem si teď prostudoval nějaké deepcopy myslím ani potřebovat nebude. Možná pro LU rozklad možná ano.

Chápu Gaussovu eliminaci dobře, když budu tvrdit, že jsou po LU rozkladu mám dvě trojůhelníkové matice a nejdůležitější na tom je to, že když počítám determinant trojůhelníkové matice, tak mi stačí počítat pouze jednu diagonálu, protože ostatní výpočty mají nejméně jeden z činitelů 0?

[Martin Tůma]

Chápu Gaussovu eliminaci dobře

Ne tak docela. LU dekompozice není částí Gaussovy eliminace. Zjednodušeně řečeno je LU dekompozice modifikovanou Gaussovou eliminací.

LU dekompozice rozkládá původní matici na součin dou trojúhelníkových matic (jedna dolní a jedna horní), Gausova eliminace matici pouze upravuje na horní trojúhelníkový tvar. K efektivnímu výpočtu determinantu se nicméně dá použít jak LU dekompozice tak Gaussova eliminace.

Jinak "finta" je skutečně v tom, že jediný nenulový součin z výpočtu determinantu je u trojúhelníkové matice právě ten na diagonále.

[Jary]

Za tři roky jsem zapoměl, co je to vlastně GE. Zatraceně.

Myslím, že už si na to začínám vzpomínat. Každopádně dík za vysvětlení.