NumPy не встает :(

[perasperaadastra]

Пытаюсь на маке установить NumPy через brew. Процесс начинается хорошо, но потом процессор загружается по полной программе на стадии "make bootstrap" и дальше я жду, жду, жду... При этом тот же numpy ставится через PyCharm за несколько секунд. В общем, homebrew явно работает не так...

[flip_flop]

Попробуйте ещё через анаконду, для разнообразия. А лучше всего собирать по старинке, из исходного кода

[perasperaadastra]

По старинке я точно не осилю... Даже homebrew для меня уже за гранью понимания. Сейчас почистил кашу и снова пробую поставить через brew. Вентилятор гудит, что-то делается, но результата пока не видно. Если бы была возможность хоть посмотреть, какие действия в данный момент совершаются...

[perasperaadastra]

Ооо! Гора родила мышь! NumPy встал! После окончания установки, выдалось совершенно издевательское сообщение:
/usr/local/Cellar/numpy/1.10.1: 1393 files, 29M, built in 2.0 minutes

Две минуты? Не знаю, по какому стандарту это две минуты, а по моим часам прошло явно больше 25.

[toly_molly]

Если бы была возможность хоть посмотреть, какие действия в данный момент совершаются...

Это как-раз легко сделать

Code: Select all

brew install homebrew/python/numpy --verbose

учитывая, что происходит полная сборка самой свежей версии numpy (включая компиляцию), то это действительно может занять довольно продолжительное время.
С --verbose флагом вы можете видеть все что происхоит (каждый файл который компилируется), так что может это привнесет утешительный эффект.

[flip_flop]

Судя по описанию процесса, HomeBrew именно и компилирует из пакетов исходников, включая все необходимые зависимости. Наберитесь терпения.

[flip_flop]

Две минуты? Не знаю, по какому стандарту это две минуты, а по моим часам прошло явно больше 25.

2 минуты на собственно NumPy а 23 минуты на пре-реквизиты.

[toly_molly]

Судя по описанию процесса, HomeBrew именно и компилирует из пакетов исходников, включая все необходимые зависимости. Наберитесь терпения.

Вот именно. Там идет полная сборка из исходников, вот поэтому я "--verbose" флаг и посоветовал, он по крайне мере будет весело на экране печатать компиляцию каждого файла исходников, всякие там "checking...", "compiling..." и прочее. Вообщем для нетерпеливых будет показывать, что компьютер бодро чего-то делает.

[perasperaadastra]

Оказалось, что паки встают правильно только с ключом --with-python3
Вроде я убил все, что ниже версии 3.5, но, видать, где-то есть какие-то остатки, раз без ключа не работает...

[perasperaadastra]

toly_molly и flip_flop: спасибо за советы и объяснения!

[flip_flop]

С --verbose флагом вы можете видеть все что происхоит (каждый файл который компилируется), так что может это привнесет утешительный эффект.

"Во многой мудрости много печали; и кто умножает познания, умножает скорбь."

Да и поздно теперь - встал, так встал. Разве что снести всё и переустановить с verbose. Но это может быть опасно для мозга

[flip_flop]

Оказалось, что паки встают правильно только с ключом --with-python3
Вроде я убил все, что ниже версии 3.5, но, видать, где-то есть какие-то остатки, раз без ключа не работает...

Python 3 & Python 2 это вообще дикое лютое веселье многие годы. Самое отвратительное в серпентарии, как мне кажется

[perasperaadastra]

С --verbose флагом вы можете видеть все что происхоит (каждый файл который компилируется), так что может это привнесет утешительный эффект.

"Во многой мудрости много печали; и кто умножает познания, умножает скорбь."

Да и поздно теперь - встал, так встал. Разве что снести всё и переустановить с verbose. Но это может быть опасно для мозга

Еще одной установки я не вынесу Это была пятая по счету. В процессе чистки после предыдущих неудачных попыток я, похоже, едва не убил систему одной неудачной командой... к счастью, ко всем системным компонентам в доступе было отказано даже с sudo.

[perasperaadastra]

Я крайне разочарован производительностью NumPy (а может быть, и самого Питона).

У меня был написан код на Processing, который симулировал столкновения двухмерных частиц. Самая ресурсоемкая часть — проверка столкновений (для всех уникальных пар частиц) и расчет их новых скоростей (операции над векторами с использованием библиотеки PVector). Получилось минимум 100 циклов в секунду (на самом деле, больше) для 400 частиц, и минимум 30 циклов в секунду для 1000 частиц

Я перевел код в Питон, полностью сохранив структуру, только адаптировал векторные операции под NumPy (через numpy.array). Результат — 2 цикла в секунду для 400 частиц! Ну как так, а?

PS Возможно, проблема еще в использовании словаря для того, чтобы обращаться к экземплярам класса через индекс (particle).

[toly_molly]

NumPy очень популярна, используется очень много где, поэтому сомневаюсь, что у вас какие-то особенные проблемы. Вполне возможно, что вы NumPy используете неоптимально, что в принципе не удивительно, если вы новичек в работе с ней. Попробуйте начать с различных статей, связанных с оптимизацией. Для затравки:
http://ipython-books.github.io/featured-01/" onclick="window.open(this.href);return false;

[flip_flop]

Python (CPython) сам по себе весьма медленный. NumPy всего лишь "обёртка" над высокопроизводительными библиотеками (BLAS, LAPACK, etc.). Теоретически, например, NumPy на основе MKL (нестандартный NumPy) должен быть примерно эквивалентен Matlab, который построен на тех же библиотеках MKL. Практически возможны нюансы из за разного использования параметров по умолчанию.

Ну и всегда есть возможность угробить производительность накладными расходами, хоть в Матлабе хоть в НумПае Есть статьи/книги/форумы на эти темы.

P.S. Прервался на поздний завтрак, а меня уже опередили Статья хорошая, кстати, применима и к Матлабу, например.

[Физик-Лирик]

Я крайне разочарован производительностью NumPy (а может быть, и самого Питона).

У меня был написан код на Processing, который симулировал столкновения двухмерных частиц. Самая ресурсоемкая часть — проверка столкновений (для всех уникальных пар частиц) и расчет их новых скоростей (операции над векторами с использованием библиотеки PVector). Получилось минимум 100 циклов в секунду (на самом деле, больше) для 400 частиц, и минимум 30 циклов в секунду для 1000 частиц

Я перевел код в Питон, полностью сохранив структуру, только адаптировал векторные операции под NumPy (через numpy.array). Результат — 2 цикла в секунду для 400 частиц! Ну как так, а?

PS Возможно, проблема еще в использовании словаря для того, чтобы обращаться к экземплярам класса через индекс (particle).

По поводу невстающей библиотеки, к сожалению, ничего сказать не могу. А не могли бы Вы добавить некоторые подробности по поводу процесса, который моделируется. Вы моделируете столкновение частиц как твердых шариков или вводите "поперечное сечение"?
Это Монте-Карловская задача, уравнение Больцмана или нечто другое?

[flip_flop]

А не могли бы Вы добавить некоторые подробности по поводу процесса, который моделируется. Вы моделируете столкновение частиц как твердых шариков или вводите "поперечное сечение"?
Это Монте-Карловская задача, уравнение Больцмана или нечто другое?

Хе-хе, старый добрый метод частиц. Вспоминается (30 лет тому назад) решение уравнения Пуассона в полупроводниках. На Фортране Метод до сих пор популярен, например для CFD, возникают новые интерпретации, типа Латтиса-Больцмана. Но для таких задач нужен Фортран, С, или, в крайнем случае, С++. А тут какой то Processing, пора заслушать начальника транспортного цеха

[perasperaadastra]

Пуассоны-Латтисы, говорите? Я и таких слов-то не знаю! У меня все просто: в свободное время я захотел познакомиться с Питоном поближе. Для ознакомления я выбрал ту же задачу, которую какое-то время назад использовал для ознакомления с Processing, а именно: смоделировать эластичное столкновение двухмерных частиц (одинаковых по размеру и массе) с помощью ньютоновских уравнений для сохранения импульса и кинэргии.

С точки зрения кода, все примитивно: есть класс Particles для частиц (id, цвет, векторные кординаты, векторная скорость). Все частицы инициализируются так:

Code: Select all

particle = {}
for i in range (0, quantity):
     particle[i] = Particles(i, (цвет), np.array(случайные координаты), np.array(случайные скорости))
     print("initialized", i)

В классе есть методы move и collision.

Метод move вызывается для каждой particle (через for i in range (0, quantity)) и перемещает каждую частицу (а также отражает ее в случае столкновения со стенками).

Метод collision тоже вызывается для каждой частицы через for i in range (0, quantity) и сравнивает положение этой частицы со всеми остальными за исключением тех, которые уже были проверены (for k in range (self.id+1, quantity))

Code: Select all

        for k in range (self.id+1, quantity):
            if (np.power(self.pos[0] + self.vel[0] - particle[k].pos[0] - particle[k].vel[0],2) + np.power(self.pos[1] + particle[k].vel[1] - particle[k].pos[1] - particle[k].vel[1],2)) <= квадрата диаметра:

В Processing это все работает на ура и получается вот так

Screen Shot 2015-12-06 at 6.30.44 PM.png

Screen Shot 2015-12-06 at 6.31.06 PM.png

Тут вам и бомбардировка, вызывающая взрывную волну, и диффузия, и даже PV=nRT.

Вот только в процессинге все интерактивно и плавно, а в гадюке все по 2 кадра в секунду.

[Физик-Лирик]

А не могли бы Вы добавить некоторые подробности по поводу процесса, который моделируется. Вы моделируете столкновение частиц как твердых шариков или вводите "поперечное сечение"?
Это Монте-Карловская задача, уравнение Больцмана или нечто другое?

Хе-хе, старый добрый метод частиц. Вспоминается (30 лет тому назад) решение уравнения Пуассона в полупроводниках. На Фортране Метод до сих пор популярен, например для CFD, возникают новые интерпретации, типа Латтиса-Больцмана. Но для таких задач нужен Фортран, С, или, в крайнем случае, С++. А тут какой то Processing, пора заслушать начальника транспортного цеха

Метод частиц 30-летней давности ... тогда Монте-Карло. Расщепление процессов на два - свободное движение
в электрическом поле, а затем "столкновение". Уравнение Пуассона .. .это два (многомерных) транспортных диффура для
электронов и дырок и уравнение Пуассона для электрического поля? Там, конечно, своих проблем хватает, т.к.
в сильных полях или в случае их больших градиентов функциия распределения частиц не является равновесной
(максвелловской), а транспортные коффициенты являются функционалами от функции распределения ... Уравнение Больцмана -
вот наш ответ.

[Физик-Лирик]

Пуассоны-Латтисы, говорите? Я и таких слов-то не знаю!

Чистое Монте-Карло. Расщепление по двум процессам: свободное движение и столкновение. Я так понимаю, столкновение моделируется по модели твердых шаров. Впрочем, если ввести дифференциальное поперечное сечение - процессы будут почти теми же. Насколько я понимаю, рассматриваются только парные столкновения - чистый Больцмановский подход. Ради фана можно попробовать доказать эквивалентность Монте-Карло и Больцмана (в данном случае), но обычно это непростая задача.
А вообще, это прикольно - изучить новый язык используя физические модели ... Снимаю шляпу. Я, увы, на такие
подвиги уже не решаюсь.

[perasperaadastra]

Кстати, если кто-нибудь захочет поиграть с моделью частиц устанавливайте Processing и вставляйте приложенный код

https://processing.org/download/?processing

Любую частицу можно выделить мышью и двигать, а чтобы затормозить все частицы, нужно нажать и держать какую-нибудь клавишу.

Code: Select all

Particles[] particle;           //object array
int quantity;
int rad;
int radsq;
float accel;
boolean selection;               // indicates if any one object already selected

void setup() {
  stroke(0);
  pixelDensity(displayDensity());
  frameRate(600);
  int jmax = 40;                 // number of particles horizontally (40)
  int kmax = 18;                 // number of particles vertically (18)
  accel = 0.000;                 // gravity acceleration (0.005)
  rad = 10;                      // radius in pixels (8) 
  int spacing = 3;               // spacing in pixels (3)
  radsq = 4 * int(sq(rad));      // radius square * 4 
  selection = false;
  
  int w = 2*jmax*(rad+spacing) + 2*spacing; // determines recommended width
  int h = 2*kmax*(rad+spacing) + 2*spacing; // determines recommended height
  print (w,h);                 // prints recommended resolution
  
  size(1046, 474);             // video frame size (1046, 474) — use w, h printed above
  
  color c = color (255, 0, 0);
  
  quantity = jmax * kmax;    //particle quantity
  
  // initialize particle objects and change the colour halfway
  particle = new Particles[quantity];
  for (int j = 0; j < jmax; j++) {
    for (int k = 0; k < kmax; k++) {
      int i = k * jmax + j;
      if (j > (jmax/2)) {
        c = color (0, 255, 0); }
      particle[i] = new Particles(c, new PVector((j*2 + 1)*(rad + spacing) + spacing, (k*2 + 1)*(rad + spacing) + spacing), new PVector(0,0)); } }
}

void draw() {
  background(255);
  for (int i = 0; i < quantity; i++) {
    particle[i].collide(i);  }             // process collission between particles 
  
  for (int i = 0; i < quantity; i++) {
    particle[i].move();                    // move and check for collision with walls
    particle[i].display();                 // display next frame
    if ((mousePressed == true) && (selection == false)) {
      particle[i].select(); }
    if (keyPressed == true) {              // slow down particles when a key is pressed
      particle[i].slowdown(); }
  }
  energy();                                // calculate and display energy
}


void energy() {                            
// calculate and display energy
  float vsum = 0;
  float hsum = 0;
  for (int i = 0; i < quantity; i++) {               
    vsum += sq(particle[i].vel.mag());               // kinetic energy
    hsum += accel * (height - particle[i].pos.y);    // gravitational energy
  }
  
  fill(0);
  textSize(20);
  text("Ek "+int(1000*vsum/quantity), 20, 20);        // kinetic energy per 1000 particles
  text("Ep "+int(1000*hsum/quantity)+"  ("+int(hsum/quantity/height/accel*200)+" %)", 20, 40);        // gravitational energy per 1000 particles and % of maximum possible
  text("E   "+int(1000*(vsum+hsum)/quantity),20, 60); // total energy per 1000 particles
}

void mouseReleased() {
  for (int i = 0; i < quantity; i++) {
    particle[i].deselect();
    selection = false;
  }
}


class Particles { 
  color c;
  PVector pos;
  PVector vel;
  boolean selected;

  Particles(color tempC, PVector temppos, PVector tempvel) { 
    c = tempC;
    pos = temppos;
    vel = tempvel;
    selected = false;
  }

  void display() {
    fill(c);
    ellipseMode(RADIUS);
    ellipse(pos.x,pos.y,rad,rad);
  }


  void collide(int j) {
  // check for and process collisions between particles (elastic)
  for (int k = j+1; k < quantity; k++) {
    if ((sq(pos.x + vel.x - particle[k].pos.x - particle[k].vel.x) +
         sq(pos.y + vel.y - particle[k].pos.y - particle[k].vel.y) <= radsq)) {
          
        PVector P = PVector.sub(particle[k].pos, pos);
        float d=P.mag();
        P.normalize();                          // unit vector in collision direction
        d -= (2*rad);
        float n = vel.mag() + particle[k].vel.mag();
        pos.add(PVector.mult(P,d*vel.mag()/(n)));
        particle[k].pos.add(PVector.mult(P,-1*d*particle[k].vel.mag()/(n)));
        // calculate speed components along and normal to collision direction
        PVector N = new PVector(P.y,P.x*(-1));  // normal vector to collision direction
        float vp1 = vel.dot(P);                 // velocity of P1 along collision direction
        float vn1 = vel.dot(N);                 // velocity of P1 normal to collision direction
        float vp2 = particle[k].vel.dot(P);     // velocity of P2 along collision direction
        float vn2 = particle[k].vel.dot(N);     // velocity of P2 normal to collision direction
        
        float vp1_after = vp2;
        float vp2_after = vp1;
        vel=PVector.add(PVector.mult(P,vp1_after),PVector.mult(N,vn1));
        particle[k].vel=PVector.add(PVector.mult(P,vp2_after),PVector.mult(N,vn2));
      }            
    }
 }

  void move() {
  // move particle and check for collision with walls
    pos.x = pos.x + vel.x;
    pos.y = pos.y + vel.y;
  
    vel.y = vel.y + accel;
    
    // check for collision with walls
    if (((pos.x + vel.x) > (width - rad)) || ((pos.x + vel.x) < rad)) {
      vel.x = -vel.x;
    }
    if (((pos.y + vel.y) > height - rad) || ((pos.y + vel.y) < rad)) {
      vel.y = -vel.y;
    }
    
    if (selected == true) {
      pos.x = mouseX;
      pos.y = mouseY;
      vel.x = mouseX - pmouseX;
      vel.y = mouseY - pmouseY;
    }
  }
  
  void select() {
    if (((pos.x - 10) < mouseX) && ((pos.x + 10) > mouseX) && ((pos.y - 10) < mouseY) && ((pos.y + 10) > mouseY) && (mousePressed == true)) {
      selected = true;
      selection = true;
    }
  }

  void slowdown() {
  //slow down particles when a key is pressed  
    vel.x *= 0.98;
    vel.y *= 0.98;
  }

  void deselect() {
    selected = false;
  }
}

[Физик-Лирик]

Спасибо за код. Я его поперек посмотрел. Идеи понятны. Несколько вопросов.
1) Каковы размеры частиц? Вы учитываете гравитационную энергию. Там что, бильярдные шары сталкиваются?
2) Во время фазы свободного движения, как Вы учитываете, столкнуться частицы или нет? Я правильно понимаю, что
свободное движение и столкновения отдельно рассматриваются? Или Вы берете время до ближайшего столкновения двух частиц,
двигаете все частицы в течение этого времени, затем парное столкновение, свободное движение, парное столкновение
и т.п.?

[perasperaadastra]

1) Сначала я планировал сделать модель идеального газа, но потом отказался от этой затеи, так как в реальном времени модель может обрабатывать очень ограниченное число частиц и с очень низкой температурой. Так что получились бильярдные шары или песчинки в зависимости от масштаба. Правда, в модели нет вращения и каких либо потерь на трение, а материал идеально упругий.

Гравитацию я добавил для того, чтобы посмотреть эффект температуры на распределение плотности по высоте. К сожалению, работает гравитация неточно (среда сама собой медленно разогревается), поэтому я ее обычно обнуляю.

2) Столкновения между частицами рассчитываются перед тем, как совершить движение. Там просто два вложенных цикла, которые проверяют каждую уникальную пару частиц на предмет наложения в следующей итерации движения. Если выясняется, что будет наложение, считается, что столкновение уже произошло и высчитываются новые скорости. Если скорости такие, что на каждой итерации частица движется не более, чем на один пиксель, то погрешность терпимая. Но если систему разогреть выше таких скоростей, то точность страдает.

После расчета столкновений для всех уникальных пар делается следующая итерация движения (для всех частиц), и снова повторяется расчет столкновений всех уникальных пар.

Эквивалентный алгоритм такой

Code: Select all

while true {
	for (j = 0; j <  quantity; j++ {
    		for (k = j+1; k < quantity; k++) { 
          		if (условие столкновения) {
          			... расчет новых скоростей } } }
          			
 // условие столкновения [(Rjx+Vjx) - (Rkx+Vkx)] + [(Rjy+Vjy) - (Rky - Vky)]^2) <= четырех радиусов в квадрате, 
 // где Rjx — x компонент радиус-вектора частицы j, а Vky — y компонент вектора скорости частицы k.

	for (j = 0; j <  quantity; j++ {
   	 	... итерация движения Rj = Rj + Vj
    		if (условие столкновения со вертикальными стенками {
   		 	 ... расчет новых скоростей } 
    		if (условие столкновения со горизонтальными стенками {
    	 		... расчет новых скоростей } }
 
	 ... вычисление энергий и вывод на экран }

[flip_flop]

Метод частиц 30-летней давности ... тогда Монте-Карло. Расщепление процессов на два - свободное движение
в электрическом поле, а затем "столкновение". Уравнение Пуассона .. .это два (многомерных) транспортных диффура для
электронов и дырок и уравнение Пуассона для электрического поля? Там, конечно, своих проблем хватает, т.к.
в сильных полях или в случае их больших градиентов функциия распределения частиц не является равновесной
(максвелловской), а транспортные коффициенты являются функционалами от функции распределения ... Уравнение Больцмана -
вот наш ответ.

Уравнение Пуассона - это образно, оно (для электрического поля) было первым в фундаментальной системе уравнений полупроводника, были ещё уравнения непрерывности для дырок и электронов, диффузия и дрейф для плотности токов , генерация-рекомбинация, распределение легирования и прочее. До К.У.Б в чистом виде редко доходили, скорее использовали его выводы, типа проводимостей или распределения носителей. Впрочем, я там стоял рядом, "не знал, не знал, и забыл". Просто в основном использовались методы конечных элементов (иногда конечных разностей), методы частиц (ну да, типа Монте Карло) выглядили иногда интересно и свежо.

Пераспераадастра, мне кажется, что в Вашей постановке доля NumPy в общем объёме вычислений мала, поэтому его влияние на конечную скорость тоже незначительно, каким бы быстрым NumPy сам по себе не был, а суммарная производительность определяется скорее чистым Питоном, который совсем не шустрый. Есть варианты более быстрого PyPy, но я не пробовал, не знаю как он дружит с NumPy.