From 46ded45310a594e3f23a0dac24fb47310f05ae4c Mon Sep 17 00:00:00 2001 From: Marcus Bannerman Date: Thu, 5 Jun 2025 15:13:00 +0100 Subject: [PATCH 01/19] Added pairwise capture state data tracking --- .../dynamo/outputplugins/collMatrix.cpp | 59 ++++++++++++++++++- .../dynamo/outputplugins/collMatrix.hpp | 28 ++++++--- 2 files changed, 76 insertions(+), 11 deletions(-) diff --git a/src/dynamo/dynamo/outputplugins/collMatrix.cpp b/src/dynamo/dynamo/outputplugins/collMatrix.cpp index f36497c4b..72f0dd1dd 100644 --- a/src/dynamo/dynamo/outputplugins/collMatrix.cpp +++ b/src/dynamo/dynamo/outputplugins/collMatrix.cpp @@ -76,8 +76,22 @@ void OPCollMatrix::eventUpdate(const Event &event, const NEventData &SDat) { auto &cs1 = _currentCaptureState[ck1]; auto &cs2 = _currentCaptureState[ck2]; - auto ek = EventKey(ck, pData.getType()); + + PairEventCaptureStateKey pecskey(ek, std::min(cs1._state, cs2._state), std::max(cs1._state, cs2._state)); + auto pit = _pairCaptureCounters.insert(decltype(_pairCaptureCounters)::value_type( + pecskey, PairEventCaptureStateData(Sim->lastRunMFT * 0.01))); + auto &pecs = pit.first->second; + if (pecs.last_event_time != 0) { + pecs.MFT.addVal(Sim->systemTime - pecs.last_event_time); + } + pecs.rijdotvij.addVal(pData.rvdot); + pecs.rijdotdP.addVal(pData.rij * pData.impulse); + pecs.vi2.addVal(pData.particle1_.getOldVel().nrm2()); + pecs.vi2.addVal(pData.particle2_.getOldVel().nrm2()); + // Now we update the last event time + pecs.last_event_time = Sim->systemTime; + auto cek1 = EventCaptureStateKey(ek, cs1._state); auto cek2 = EventCaptureStateKey(ek, cs2._state); @@ -271,8 +285,47 @@ void OPCollMatrix::output(magnet::xml::XmlStream &XML) { XML << magnet::xml::endtag("Count"); } - XML << magnet::xml::endtag("CaptureCounters") - << magnet::xml::tag("CaptureStateHistogram"); + XML << magnet::xml::endtag("CaptureCounters"); + + XML << magnet::xml::tag("PairCaptureCounters"); + for (const auto &val : _pairCaptureCounters) { + auto cek = val.first; + auto ek = std::get<0>(cek); + auto class_key = ek.first; + auto event_type = ek.second; + auto captures1 = std::get<1>(cek); + auto captures2 = std::get<2>(cek); + auto &pecs = val.second; + + XML << magnet::xml::tag("Count") << magnet::xml::attr("Name") + << getEventSourceName(class_key, Sim) << magnet::xml::attr("Event") + << event_type << magnet::xml::attr("captures1") << captures1 + << magnet::xml::attr("captures2") << captures2; + + XML << magnet::xml::tag("MFT"); + pecs.MFT.outputHistogram(XML, 1.0 / Sim->units.unitTime()); + XML << magnet::xml::endtag("MFT"); + + XML << magnet::xml::tag("RijDotVij"); + pecs.rijdotvij.outputHistogram(XML, 1.0 / Sim->units.unitLength() / + Sim->units.unitVelocity()); + XML << magnet::xml::endtag("RijDotVij"); + + XML << magnet::xml::tag("RijDotDeltaPij"); + pecs.rijdotdP.outputHistogram(XML, 1.0 / Sim->units.unitLength() / + Sim->units.unitMomentum()); + XML << magnet::xml::endtag("RijDotDeltaPij"); + + XML << magnet::xml::tag("V2"); + pecs.vi2.outputHistogram(XML, 1.0 / Sim->units.unitVelocity() / + Sim->units.unitVelocity()); + XML << magnet::xml::endtag("V2"); + + XML << magnet::xml::endtag("Count"); + } + XML << magnet::xml::endtag("PairCaptureCounters"); + + XML << magnet::xml::tag("CaptureStateHistogram"); // Before we output the histogram we need to bring everything up to date for (auto &p : _currentCaptureState) { diff --git a/src/dynamo/dynamo/outputplugins/collMatrix.hpp b/src/dynamo/dynamo/outputplugins/collMatrix.hpp index 45d5ee6a8..c967e4529 100644 --- a/src/dynamo/dynamo/outputplugins/collMatrix.hpp +++ b/src/dynamo/dynamo/outputplugins/collMatrix.hpp @@ -51,18 +51,17 @@ class OPCollMatrix : public OutputPlugin { unsigned long totalCount; - // We create a key for events based on the interaction/system/global/local ID - // and type (EventSourceKey) and EventType + // EventKet is a pair of EventSourceKey and EEventType + // It describes the type of event and its source //! A key for two events - // Used to track the previous and current event for some data typedef std::pair InterEventKey; + //! Counters for properties between two different events std::map counters; // First we track how many times a particle has been captured - typedef std::pair - TotalCaptureStateKey; // Interaction ID and particle ID + typedef std::pair TotalCaptureStateKey; // Interaction ID and particle ID struct CaptureStateData { CaptureStateData(double binWidth = 1.0) {} double _last_update = 0; @@ -76,9 +75,8 @@ class OPCollMatrix : public OutputPlugin { magnet::math::HistogramWeighted<> _captureStateHistogram; // Here we're tracking collision statistics depending on the Event Type/Source - // and pair capture state + // and capture state of each particle individually typedef std::pair EventCaptureStateKey; - struct EventCaptureStateData { EventCaptureStateData(double binWidth) : MFT(binWidth), rijdotvij(0.01), rijdotdP(0.01), vi2(0.01) {} @@ -89,9 +87,23 @@ class OPCollMatrix : public OutputPlugin { magnet::math::Histogram<> vi2; magnet::math::Histogram<> _particle_MFT; }; - std::map _captureCounters; + typedef std::tuple PairEventCaptureStateKey; + struct PairEventCaptureStateData { + PairEventCaptureStateData(double binWidth) + : MFT(binWidth), rijdotvij(0.01), rijdotdP(0.01), vi2(0.01) {} + double last_event_time = 0; + magnet::math::Histogram<> MFT; + magnet::math::Histogram<> rijdotvij; + magnet::math::Histogram<> rijdotdP; + magnet::math::Histogram<> vi2; + magnet::math::Histogram<> _particle_MFT; + }; + std::map _pairCaptureCounters; + + + typedef std::pair MFTKey; std::map> _fullMFT; From 9ebeee86d6800e278cb6820476412e1a5b5fd0da Mon Sep 17 00:00:00 2001 From: Marcus Bannerman Date: Tue, 29 Apr 2025 13:34:53 +0100 Subject: [PATCH 02/19] Creating a Chung-Lu plotting script --- scripts/SW_eos.py | 4 +- scripts/plotter.ipynb | 1019 ++++++++++++++++++++++++++++++++++++++ src/pydynamo/__init__.py | 16 +- 3 files changed, 1028 insertions(+), 11 deletions(-) create mode 100644 scripts/plotter.ipynb diff --git a/scripts/SW_eos.py b/scripts/SW_eos.py index c71ba7e87..894aa22b6 100755 --- a/scripts/SW_eos.py +++ b/scripts/SW_eos.py @@ -151,6 +151,6 @@ def setup_worker( config, # The name of the config file to genera # This creates a pandas dataframe with columns for the state variables # AND any output values. It also generates pkl files, some for # different properties. -data = mgr.fetch_data(1000) +df, state_data = mgr.fetch_data(1000) -print(data) +print(df) diff --git a/scripts/plotter.ipynb b/scripts/plotter.ipynb new file mode 100644 index 000000000..3b1a19f68 --- /dev/null +++ b/scripts/plotter.ipynb @@ -0,0 +1,1019 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 53, + "id": "bbb1036b", + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
NEventsTottTotalpcvu
InitStateLambdakTndensityN
FCC1.51.00.0140004000000023845.0359320.009358+/-0.0000090.148+/-0.004-0.13712+/-0.00019
0.104000400000001247.1170080.0319+/-0.000517+/-6-2.57+/-0.06
0.50400040000000682.734445-0.039+/-0.0043.2+/-0.6-4.648+/-0.013
1.00400040000000352.0543658.353+/-0.0070.284+/-0.010-6.9236+/-0.0008
1.30400040000000153.40526040.61+/-0.050.03097+/-0.00030-8.96942+/-0.00004
1.50.0140004000000025590.7012910.014634+/-0.0000110.0441+/-0.0013-0.09652+/-0.00014
0.104000400000002671.5543740.11652+/-0.000160.426+/-0.020-0.9247+/-0.0009
0.50400040000000655.9841670.5409+/-0.00210.407+/-0.018-3.7988+/-0.0009
1.00400040000000299.92020113.437+/-0.0100.1274+/-0.0020-6.82028+/-0.00023
1.30400040000000121.66455963.43+/-0.060.02037+/-0.00027-8.95679+/-0.00005
2.00.0140004000000025019.7893740.019776+/-0.0000160.0202+/-0.0004-0.08169+/-0.00013
0.104000400000002599.1760140.18250+/-0.000150.171+/-0.005-0.7881+/-0.0005
0.50400040000000574.3013801.3371+/-0.00330.185+/-0.004-3.6622+/-0.0006
1.00400040000000265.20732018.582+/-0.0190.0699+/-0.0020-6.77118+/-0.00028
1.30400040000000103.53887986.59+/-0.100.01361+/-0.00016-8.94862+/-0.00004
2.50.0140004000000023962.1589450.024930+/-0.0000230.01144+/-0.00015-0.07403+/-0.00011
0.104000400000002460.5163920.24633+/-0.000320.0953+/-0.0018-0.7248+/-0.0005
0.50400040000000516.0677962.1419+/-0.00240.114+/-0.004-3.58867+/-0.00035
1.00400040000000240.08835823.714+/-0.0190.0465+/-0.0015-6.74276+/-0.00026
1.3040004000000091.497844110.05+/-0.100.00965+/-0.00009-8.94287+/-0.00006
3.00.0140004000000022859.6378920.030053+/-0.0000230.00768+/-0.00009-0.06928+/-0.00005
0.104000400000002326.2656720.30937+/-0.000280.0581+/-0.0012-0.68636+/-0.00024
0.50400040000000471.6295082.9577+/-0.00320.0778+/-0.0020-3.54119+/-0.00031
1.00400040000000220.99661328.857+/-0.0270.0330+/-0.0007-6.72278+/-0.00026
1.3040004000000082.885040133.37+/-0.130.00727+/-0.00010-8.93878+/-0.00009
inf0.0140004000000048493.5669770.010216+/-0.0000100.0+/-00.0+/-0
0.104000400000004638.7246700.12396+/-0.000110.0+/-00.0+/-0
0.50400040000000814.1862161.6356+/-0.00100.0+/-00.0+/-0
1.00400040000000398.44075510.252+/-0.0100.0+/-00.0+/-0
1.30400040000000136.30320747.68+/-0.040.0+/-00.0+/-0
2.01.00.01400040000000461.9272110.00368+/-0.00021(4.1+/-1.0)e+02-6.0+/-0.4
0.10400040000000198.664502-0.0025+/-0.0010(2.0+/-0.7)e+02-13.17+/-0.29
0.50400040000000143.602432-0.646+/-0.01886+/-20-16.17+/-0.18
1.00400040000000164.077613-28.683+/-0.0072.08+/-0.15-19.9221+/-0.0016
1.30400040000000150.88231547.854+/-0.0340.001190+/-0.000019-21.001139+/-0.000010
1.50.0140004000000015156.3574710.013198+/-0.0000100.1713+/-0.0029-0.2990+/-0.0004
0.10400040000000256.3329240.0033+/-0.000435+/-15-9.76+/-0.16
0.50400040000000173.414606-0.400+/-0.02740+/-17-13.36+/-0.14
1.00400040000000125.644637-25.755+/-0.0191.65+/-0.09-18.8927+/-0.0017
1.30400040000000124.82870571.75+/-0.070.000360+/-0.000004-21.0008050+/-0.0000035
2.00.0140004000000015303.8735260.018525+/-0.0000210.0692+/-0.0016-0.24377+/-0.00024
0.10400040000000448.0477390.0400+/-0.001331+/-10-6.12+/-0.19
0.50400040000000269.417239-0.231+/-0.0178.9+/-2.9-9.67+/-0.08
1.00400040000000111.497971-18.400+/-0.0310.74+/-0.04-18.3158+/-0.0018
1.30400040000000108.68044195.57+/-0.080.0001729+/-0.0000011-21.0006807+/-0.0000027
2.50.0140004000000014762.0702270.023735+/-0.0000130.0371+/-0.0005-0.21868+/-0.00021
0.104000400000001410.9479040.13783+/-0.000300.74+/-0.05-2.204+/-0.004
0.50400040000000311.9337870.238+/-0.0040.263+/-0.019-8.3044+/-0.0013
1.00400040000000101.762522-11.65+/-0.040.448+/-0.015-18.0130+/-0.0015
1.3040004000000097.506480119.37+/-0.100.0001011+/-0.0000007-21.0006130+/-0.0000022
3.00.0140004000000014116.9469080.028891+/-0.0000250.02296+/-0.00028-0.20397+/-0.00016
0.104000400000001426.6386580.2057+/-0.00050.280+/-0.013-1.9688+/-0.0020
0.50400040000000290.4862491.0408+/-0.00320.162+/-0.010-8.2026+/-0.0006
1.0040004000000093.828259-5.50+/-0.040.310+/-0.007-17.8273+/-0.0012
1.3040004000000089.105293143.35+/-0.06(6.44+/-0.04)e-05-21.0005735+/-0.0000031
inf0.0140004000000029776.4998140.0101969+/-0.00000340.0+/-00.0+/-0
0.104000400000002944.7749580.12377+/-0.000120.0+/-00.0+/-0
0.50400040000000532.2841541.6372+/-0.00120.0+/-00.0+/-0
1.00400040000000165.81394010.247+/-0.0110.0+/-00.0+/-0
1.30400040000000155.36996947.78+/-0.040.0+/-00.0+/-0
\n", + "
" + ], + "text/plain": [ + " NEventsTot tTotal \\\n", + "InitState Lambda kT ndensity N \n", + "FCC 1.5 1.0 0.01 4000 40000000 23845.035932 \n", + " 0.10 4000 40000000 1247.117008 \n", + " 0.50 4000 40000000 682.734445 \n", + " 1.00 4000 40000000 352.054365 \n", + " 1.30 4000 40000000 153.405260 \n", + " 1.5 0.01 4000 40000000 25590.701291 \n", + " 0.10 4000 40000000 2671.554374 \n", + " 0.50 4000 40000000 655.984167 \n", + " 1.00 4000 40000000 299.920201 \n", + " 1.30 4000 40000000 121.664559 \n", + " 2.0 0.01 4000 40000000 25019.789374 \n", + " 0.10 4000 40000000 2599.176014 \n", + " 0.50 4000 40000000 574.301380 \n", + " 1.00 4000 40000000 265.207320 \n", + " 1.30 4000 40000000 103.538879 \n", + " 2.5 0.01 4000 40000000 23962.158945 \n", + " 0.10 4000 40000000 2460.516392 \n", + " 0.50 4000 40000000 516.067796 \n", + " 1.00 4000 40000000 240.088358 \n", + " 1.30 4000 40000000 91.497844 \n", + " 3.0 0.01 4000 40000000 22859.637892 \n", + " 0.10 4000 40000000 2326.265672 \n", + " 0.50 4000 40000000 471.629508 \n", + " 1.00 4000 40000000 220.996613 \n", + " 1.30 4000 40000000 82.885040 \n", + " inf 0.01 4000 40000000 48493.566977 \n", + " 0.10 4000 40000000 4638.724670 \n", + " 0.50 4000 40000000 814.186216 \n", + " 1.00 4000 40000000 398.440755 \n", + " 1.30 4000 40000000 136.303207 \n", + " 2.0 1.0 0.01 4000 40000000 461.927211 \n", + " 0.10 4000 40000000 198.664502 \n", + " 0.50 4000 40000000 143.602432 \n", + " 1.00 4000 40000000 164.077613 \n", + " 1.30 4000 40000000 150.882315 \n", + " 1.5 0.01 4000 40000000 15156.357471 \n", + " 0.10 4000 40000000 256.332924 \n", + " 0.50 4000 40000000 173.414606 \n", + " 1.00 4000 40000000 125.644637 \n", + " 1.30 4000 40000000 124.828705 \n", + " 2.0 0.01 4000 40000000 15303.873526 \n", + " 0.10 4000 40000000 448.047739 \n", + " 0.50 4000 40000000 269.417239 \n", + " 1.00 4000 40000000 111.497971 \n", + " 1.30 4000 40000000 108.680441 \n", + " 2.5 0.01 4000 40000000 14762.070227 \n", + " 0.10 4000 40000000 1410.947904 \n", + " 0.50 4000 40000000 311.933787 \n", + " 1.00 4000 40000000 101.762522 \n", + " 1.30 4000 40000000 97.506480 \n", + " 3.0 0.01 4000 40000000 14116.946908 \n", + " 0.10 4000 40000000 1426.638658 \n", + " 0.50 4000 40000000 290.486249 \n", + " 1.00 4000 40000000 93.828259 \n", + " 1.30 4000 40000000 89.105293 \n", + " inf 0.01 4000 40000000 29776.499814 \n", + " 0.10 4000 40000000 2944.774958 \n", + " 0.50 4000 40000000 532.284154 \n", + " 1.00 4000 40000000 165.813940 \n", + " 1.30 4000 40000000 155.369969 \n", + "\n", + " p \\\n", + "InitState Lambda kT ndensity N \n", + "FCC 1.5 1.0 0.01 4000 0.009358+/-0.000009 \n", + " 0.10 4000 0.0319+/-0.0005 \n", + " 0.50 4000 -0.039+/-0.004 \n", + " 1.00 4000 8.353+/-0.007 \n", + " 1.30 4000 40.61+/-0.05 \n", + " 1.5 0.01 4000 0.014634+/-0.000011 \n", + " 0.10 4000 0.11652+/-0.00016 \n", + " 0.50 4000 0.5409+/-0.0021 \n", + " 1.00 4000 13.437+/-0.010 \n", + " 1.30 4000 63.43+/-0.06 \n", + " 2.0 0.01 4000 0.019776+/-0.000016 \n", + " 0.10 4000 0.18250+/-0.00015 \n", + " 0.50 4000 1.3371+/-0.0033 \n", + " 1.00 4000 18.582+/-0.019 \n", + " 1.30 4000 86.59+/-0.10 \n", + " 2.5 0.01 4000 0.024930+/-0.000023 \n", + " 0.10 4000 0.24633+/-0.00032 \n", + " 0.50 4000 2.1419+/-0.0024 \n", + " 1.00 4000 23.714+/-0.019 \n", + " 1.30 4000 110.05+/-0.10 \n", + " 3.0 0.01 4000 0.030053+/-0.000023 \n", + " 0.10 4000 0.30937+/-0.00028 \n", + " 0.50 4000 2.9577+/-0.0032 \n", + " 1.00 4000 28.857+/-0.027 \n", + " 1.30 4000 133.37+/-0.13 \n", + " inf 0.01 4000 0.010216+/-0.000010 \n", + " 0.10 4000 0.12396+/-0.00011 \n", + " 0.50 4000 1.6356+/-0.0010 \n", + " 1.00 4000 10.252+/-0.010 \n", + " 1.30 4000 47.68+/-0.04 \n", + " 2.0 1.0 0.01 4000 0.00368+/-0.00021 \n", + " 0.10 4000 -0.0025+/-0.0010 \n", + " 0.50 4000 -0.646+/-0.018 \n", + " 1.00 4000 -28.683+/-0.007 \n", + " 1.30 4000 47.854+/-0.034 \n", + " 1.5 0.01 4000 0.013198+/-0.000010 \n", + " 0.10 4000 0.0033+/-0.0004 \n", + " 0.50 4000 -0.400+/-0.027 \n", + " 1.00 4000 -25.755+/-0.019 \n", + " 1.30 4000 71.75+/-0.07 \n", + " 2.0 0.01 4000 0.018525+/-0.000021 \n", + " 0.10 4000 0.0400+/-0.0013 \n", + " 0.50 4000 -0.231+/-0.017 \n", + " 1.00 4000 -18.400+/-0.031 \n", + " 1.30 4000 95.57+/-0.08 \n", + " 2.5 0.01 4000 0.023735+/-0.000013 \n", + " 0.10 4000 0.13783+/-0.00030 \n", + " 0.50 4000 0.238+/-0.004 \n", + " 1.00 4000 -11.65+/-0.04 \n", + " 1.30 4000 119.37+/-0.10 \n", + " 3.0 0.01 4000 0.028891+/-0.000025 \n", + " 0.10 4000 0.2057+/-0.0005 \n", + " 0.50 4000 1.0408+/-0.0032 \n", + " 1.00 4000 -5.50+/-0.04 \n", + " 1.30 4000 143.35+/-0.06 \n", + " inf 0.01 4000 0.0101969+/-0.0000034 \n", + " 0.10 4000 0.12377+/-0.00012 \n", + " 0.50 4000 1.6372+/-0.0012 \n", + " 1.00 4000 10.247+/-0.011 \n", + " 1.30 4000 47.78+/-0.04 \n", + "\n", + " cv \\\n", + "InitState Lambda kT ndensity N \n", + "FCC 1.5 1.0 0.01 4000 0.148+/-0.004 \n", + " 0.10 4000 17+/-6 \n", + " 0.50 4000 3.2+/-0.6 \n", + " 1.00 4000 0.284+/-0.010 \n", + " 1.30 4000 0.03097+/-0.00030 \n", + " 1.5 0.01 4000 0.0441+/-0.0013 \n", + " 0.10 4000 0.426+/-0.020 \n", + " 0.50 4000 0.407+/-0.018 \n", + " 1.00 4000 0.1274+/-0.0020 \n", + " 1.30 4000 0.02037+/-0.00027 \n", + " 2.0 0.01 4000 0.0202+/-0.0004 \n", + " 0.10 4000 0.171+/-0.005 \n", + " 0.50 4000 0.185+/-0.004 \n", + " 1.00 4000 0.0699+/-0.0020 \n", + " 1.30 4000 0.01361+/-0.00016 \n", + " 2.5 0.01 4000 0.01144+/-0.00015 \n", + " 0.10 4000 0.0953+/-0.0018 \n", + " 0.50 4000 0.114+/-0.004 \n", + " 1.00 4000 0.0465+/-0.0015 \n", + " 1.30 4000 0.00965+/-0.00009 \n", + " 3.0 0.01 4000 0.00768+/-0.00009 \n", + " 0.10 4000 0.0581+/-0.0012 \n", + " 0.50 4000 0.0778+/-0.0020 \n", + " 1.00 4000 0.0330+/-0.0007 \n", + " 1.30 4000 0.00727+/-0.00010 \n", + " inf 0.01 4000 0.0+/-0 \n", + " 0.10 4000 0.0+/-0 \n", + " 0.50 4000 0.0+/-0 \n", + " 1.00 4000 0.0+/-0 \n", + " 1.30 4000 0.0+/-0 \n", + " 2.0 1.0 0.01 4000 (4.1+/-1.0)e+02 \n", + " 0.10 4000 (2.0+/-0.7)e+02 \n", + " 0.50 4000 86+/-20 \n", + " 1.00 4000 2.08+/-0.15 \n", + " 1.30 4000 0.001190+/-0.000019 \n", + " 1.5 0.01 4000 0.1713+/-0.0029 \n", + " 0.10 4000 35+/-15 \n", + " 0.50 4000 40+/-17 \n", + " 1.00 4000 1.65+/-0.09 \n", + " 1.30 4000 0.000360+/-0.000004 \n", + " 2.0 0.01 4000 0.0692+/-0.0016 \n", + " 0.10 4000 31+/-10 \n", + " 0.50 4000 8.9+/-2.9 \n", + " 1.00 4000 0.74+/-0.04 \n", + " 1.30 4000 0.0001729+/-0.0000011 \n", + " 2.5 0.01 4000 0.0371+/-0.0005 \n", + " 0.10 4000 0.74+/-0.05 \n", + " 0.50 4000 0.263+/-0.019 \n", + " 1.00 4000 0.448+/-0.015 \n", + " 1.30 4000 0.0001011+/-0.0000007 \n", + " 3.0 0.01 4000 0.02296+/-0.00028 \n", + " 0.10 4000 0.280+/-0.013 \n", + " 0.50 4000 0.162+/-0.010 \n", + " 1.00 4000 0.310+/-0.007 \n", + " 1.30 4000 (6.44+/-0.04)e-05 \n", + " inf 0.01 4000 0.0+/-0 \n", + " 0.10 4000 0.0+/-0 \n", + " 0.50 4000 0.0+/-0 \n", + " 1.00 4000 0.0+/-0 \n", + " 1.30 4000 0.0+/-0 \n", + "\n", + " u \n", + "InitState Lambda kT ndensity N \n", + "FCC 1.5 1.0 0.01 4000 -0.13712+/-0.00019 \n", + " 0.10 4000 -2.57+/-0.06 \n", + " 0.50 4000 -4.648+/-0.013 \n", + " 1.00 4000 -6.9236+/-0.0008 \n", + " 1.30 4000 -8.96942+/-0.00004 \n", + " 1.5 0.01 4000 -0.09652+/-0.00014 \n", + " 0.10 4000 -0.9247+/-0.0009 \n", + " 0.50 4000 -3.7988+/-0.0009 \n", + " 1.00 4000 -6.82028+/-0.00023 \n", + " 1.30 4000 -8.95679+/-0.00005 \n", + " 2.0 0.01 4000 -0.08169+/-0.00013 \n", + " 0.10 4000 -0.7881+/-0.0005 \n", + " 0.50 4000 -3.6622+/-0.0006 \n", + " 1.00 4000 -6.77118+/-0.00028 \n", + " 1.30 4000 -8.94862+/-0.00004 \n", + " 2.5 0.01 4000 -0.07403+/-0.00011 \n", + " 0.10 4000 -0.7248+/-0.0005 \n", + " 0.50 4000 -3.58867+/-0.00035 \n", + " 1.00 4000 -6.74276+/-0.00026 \n", + " 1.30 4000 -8.94287+/-0.00006 \n", + " 3.0 0.01 4000 -0.06928+/-0.00005 \n", + " 0.10 4000 -0.68636+/-0.00024 \n", + " 0.50 4000 -3.54119+/-0.00031 \n", + " 1.00 4000 -6.72278+/-0.00026 \n", + " 1.30 4000 -8.93878+/-0.00009 \n", + " inf 0.01 4000 0.0+/-0 \n", + " 0.10 4000 0.0+/-0 \n", + " 0.50 4000 0.0+/-0 \n", + " 1.00 4000 0.0+/-0 \n", + " 1.30 4000 0.0+/-0 \n", + " 2.0 1.0 0.01 4000 -6.0+/-0.4 \n", + " 0.10 4000 -13.17+/-0.29 \n", + " 0.50 4000 -16.17+/-0.18 \n", + " 1.00 4000 -19.9221+/-0.0016 \n", + " 1.30 4000 -21.001139+/-0.000010 \n", + " 1.5 0.01 4000 -0.2990+/-0.0004 \n", + " 0.10 4000 -9.76+/-0.16 \n", + " 0.50 4000 -13.36+/-0.14 \n", + " 1.00 4000 -18.8927+/-0.0017 \n", + " 1.30 4000 -21.0008050+/-0.0000035 \n", + " 2.0 0.01 4000 -0.24377+/-0.00024 \n", + " 0.10 4000 -6.12+/-0.19 \n", + " 0.50 4000 -9.67+/-0.08 \n", + " 1.00 4000 -18.3158+/-0.0018 \n", + " 1.30 4000 -21.0006807+/-0.0000027 \n", + " 2.5 0.01 4000 -0.21868+/-0.00021 \n", + " 0.10 4000 -2.204+/-0.004 \n", + " 0.50 4000 -8.3044+/-0.0013 \n", + " 1.00 4000 -18.0130+/-0.0015 \n", + " 1.30 4000 -21.0006130+/-0.0000022 \n", + " 3.0 0.01 4000 -0.20397+/-0.00016 \n", + " 0.10 4000 -1.9688+/-0.0020 \n", + " 0.50 4000 -8.2026+/-0.0006 \n", + " 1.00 4000 -17.8273+/-0.0012 \n", + " 1.30 4000 -21.0005735+/-0.0000031 \n", + " inf 0.01 4000 0.0+/-0 \n", + " 0.10 4000 0.0+/-0 \n", + " 0.50 4000 0.0+/-0 \n", + " 1.00 4000 0.0+/-0 \n", + " 1.30 4000 0.0+/-0 " + ] + }, + "execution_count": 53, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "import pandas\n", + "import pickle\n", + "import ipecharts\n", + "import pydynamo\n", + "import uncertainties\n", + "import collections\n", + "import matplotlib\n", + "import matplotlib.pyplot as plt\n", + "\n", + "\n", + "base = \"SW_eos\"\n", + "# Load the files from the fetch_data\n", + "state_data = pickle.load(open(f\"{base}.raw_data.pkl\", \"rb\"))\n", + "state_vars = set([j[0] for k in state_data.keys() for j in k])\n", + "\n", + "# Turn state variables into a multi-index, with a particular preferred order\n", + "preferred_state_order = (\"Lambda\", \"kT\", \"ndensity\", \"N\")\n", + "state_vars = sorted(state_vars, key=lambda x: preferred_state_order.index(x) if x in preferred_state_order else -1)\n", + "orig_df = pickle.load(open(f\"{base}.df.pkl\", \"rb\"))\n", + "orig_df.set_index(list(state_vars), inplace=True)\n", + "\n", + "# Filter out columns that do not change (boring state variables i.e.)\n", + "#orig_df = orig_df.loc[:, orig_df.nunique() > 1]\n", + "\n", + "\n", + "# Create a dataframe of scalars\n", + "scalar_df = orig_df.copy()\n", + "for col in scalar_df.columns:\n", + " if scalar_df[col].dtype == \"object\":\n", + " if isinstance(scalar_df[col].iloc[0], collections.abc.Iterable):\n", + " del scalar_df[col]\n", + "scalar_df" + ] + }, + { + "cell_type": "code", + "execution_count": 55, + "id": "3c5f65e0", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "InitState Lambda kT ndensity N \n", + "FCC 1.5 1.0 0.01 4000 [(1, 1), (2, 2), (1, 2), (1, 3), (4, 5), (3, 5...\n", + " 0.10 4000 [(8, 10), (10, 10), (6, 10), (10, 11), (7, 10)...\n", + " 0.50 4000 [(8, 9), (8, 10), (8, 8), (8, 11), (8, 12), (4...\n", + " 1.00 4000 [(13, 16), (13, 14), (13, 13), (12, 13), (14, ...\n", + " 1.30 4000 [(18, 18), (17, 18), (17, 17), (16, 18), (16, ...\n", + " 1.5 0.01 4000 [(1, 1), (2, 3), (1, 3), (1, 2), (2, 2), (2, 4...\n", + " 0.10 4000 [(1, 2), (1, 1), (2, 4), (2, 3), (3, 3), (3, 5...\n", + " 0.50 4000 [(6, 7), (6, 8), (5, 6), (6, 9), (7, 8), (8, 1...\n", + " 1.00 4000 [(12, 14), (13, 14), (14, 14), (14, 15), (13, ...\n", + " 1.30 4000 [(18, 18), (17, 18), (16, 18), (17, 17), (16, ...\n", + " 2.0 0.01 4000 [(1, 1), (2, 2), (1, 2), (1, 3), (2, 3), (3, 3)]\n", + " 0.10 4000 [(1, 3), (2, 4), (2, 3), (2, 2), (4, 6), (4, 4...\n", + " 0.50 4000 [(7, 8), (8, 10), (8, 9), (8, 8), (6, 8), (6, ...\n", + " 1.00 4000 [(14, 14), (14, 15), (13, 14), (12, 14), (13, ...\n", + " 1.30 4000 [(17, 18), (18, 18), (16, 18), (17, 17), (16, ...\n", + " 2.5 0.01 4000 [(1, 1), (1, 2), (2, 2), (1, 3), (2, 3), (3, 3)]\n", + " 0.10 4000 [(1, 2), (1, 3), (2, 3), (1, 1), (3, 4), (4, 4...\n", + " 0.50 4000 [(5, 5), (5, 7), (5, 8), (8, 8), (8, 10), (7, ...\n", + " 1.00 4000 [(14, 14), (13, 14), (12, 14), (14, 15), (12, ...\n", + " 1.30 4000 [(18, 18), (17, 18), (16, 18), (17, 17), (16, ...\n", + " 3.0 0.01 4000 [(1, 2), (1, 1), (2, 2), (1, 3), (2, 3), (3, 3)]\n", + " 0.10 4000 [(2, 3), (1, 1), (1, 2), (2, 4), (4, 4), (1, 3...\n", + " 0.50 4000 [(9, 10), (8, 9), (9, 11), (7, 9), (6, 9), (9,...\n", + " 1.00 4000 [(12, 13), (12, 14), (12, 15), (13, 14), (14, ...\n", + " 1.30 4000 [(18, 18), (17, 18), (17, 17), (16, 18), (16, ...\n", + " inf 0.01 4000 [(1, 1), (1, 2), (2, 2), (2, 3), (1, 3), (3, 3)]\n", + " 0.10 4000 [(2, 2), (1, 2), (2, 3), (1, 3), (1, 4), (1, 1...\n", + " 0.50 4000 [(6, 6), (5, 6), (6, 10), (6, 8), (5, 7), (6, ...\n", + " 1.00 4000 [(13, 15), (13, 13), (13, 14), (11, 13), (12, ...\n", + " 1.30 4000 [(18, 18), (17, 18), (16, 18), (17, 17), (16, ...\n", + " 2.0 1.0 0.01 4000 [(21, 32), (21, 48), (21, 41), (21, 34), (21, ...\n", + " 0.10 4000 [(38, 38), (37, 38), (22, 38), (24, 38), (19, ...\n", + " 0.50 4000 [(34, 45), (45, 45), (43, 45), (36, 45), (22, ...\n", + " 1.00 4000 [(39, 40), (40, 41), (40, 42), (40, 40), (38, ...\n", + " 1.30 4000 [(42, 42), (42, 43), (43, 43)]\n", + " 1.5 0.01 4000 [(1, 1), (1, 2), (2, 2), (2, 3), (4, 4), (3, 3...\n", + " 0.10 4000 [(12, 14), (12, 13), (12, 20), (12, 19), (12, ...\n", + " 0.50 4000 [(29, 31), (31, 32), (31, 31), (24, 31), (27, ...\n", + " 1.00 4000 [(38, 38), (33, 38), (36, 38), (35, 38), (38, ...\n", + " 1.30 4000 [(42, 42), (42, 43), (43, 43)]\n", + " 2.0 0.01 4000 [(1, 1), (2, 3), (1, 3), (3, 3), (2, 2), (1, 2...\n", + " 0.10 4000 [(2, 2), (22, 22), (21, 22), (22, 24), (15, 22...\n", + " 0.50 4000 [(18, 21), (17, 18), (18, 19), (18, 20), (16, ...\n", + " 1.00 4000 [(34, 41), (34, 39), (34, 37), (34, 38), (34, ...\n", + " 1.30 4000 [(42, 42), (42, 43), (43, 43)]\n", + " 2.5 0.01 4000 [(1, 2), (1, 3), (1, 1), (2, 2), (2, 4), (3, 4...\n", + " 0.10 4000 [(4, 8), (4, 6), (3, 4), (6, 11), (5, 6), (6, ...\n", + " 0.50 4000 [(18, 20), (17, 20), (20, 21), (19, 20), (20, ...\n", + " 1.00 4000 [(31, 37), (37, 38), (35, 37), (36, 37), (34, ...\n", + " 1.30 4000 [(42, 42), (42, 43), (43, 43)]\n", + " 3.0 0.01 4000 [(1, 3), (1, 1), (1, 2), (2, 3), (2, 2), (3, 3...\n", + " 0.10 4000 [(5, 8), (8, 8), (6, 8), (7, 8), (1, 3), (2, 6...\n", + " 0.50 4000 [(16, 16), (15, 16), (16, 17), (16, 18), (16, ...\n", + " 1.00 4000 [(36, 38), (35, 38), (38, 39), (38, 38), (38, ...\n", + " 1.30 4000 [(42, 42), (42, 43), (43, 43)]\n", + " inf 0.01 4000 [(1, 1), (1, 3), (2, 2), (2, 3), (1, 2), (3, 3...\n", + " 0.10 4000 [(3, 4), (3, 3), (2, 3), (5, 6), (6, 7), (3, 6...\n", + " 0.50 4000 [(13, 16), (15, 16), (16, 17), (14, 16), (16, ...\n", + " 1.00 4000 [(36, 37), (36, 36), (35, 36), (31, 36), (30, ...\n", + " 1.30 4000 [(42, 42), (42, 43), (43, 43)]\n", + "Name: ChungLu, dtype: object" + ] + }, + "execution_count": 55, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "orig_df[\"ChungLu\"]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "acb5456f", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": ".venv", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.13.2" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/src/pydynamo/__init__.py b/src/pydynamo/__init__.py index b7c9a8db9..e3fa0fc25 100755 --- a/src/pydynamo/__init__.py +++ b/src/pydynamo/__init__.py @@ -133,7 +133,7 @@ def worker(state, workdir, outputplugins, particle_equil_events, particle_run_ev except subprocess.CalledProcessError as e: raise RuntimeError('Failed while running worker, command was\n"'+str(e.cmd)+'"\nSee logfile "'+str(os.path.join(workdir, 'run.log'))+'"') -def perdir(args): +def fetch_data_worker(args): output_dir, particle_equil_events, manager = args output_dir = os.path.join(manager.workdir, output_dir) if not os.path.isdir(output_dir): @@ -518,13 +518,13 @@ def fetch_data(self, particle_equil_events, only_current_statevars = False): #We store the extracted data in a dict of dicts. The first #dict is for the state, the second for the property. - state_data = collections.defaultdict(dict) + state_data = collections.defaultdict(dict) - #So we run the per data dir operation, then reduce everything + #So we run the per data dir operation, then reduce everything. state_data = {} with alive_progress.alive_bar(n) as progress: #This is a parallel loop, returning items as they finish in arbitrary order - for result in self.imap_unordered(perdir, [(d, particle_equil_events, self) for d in output_dirs]): + for result in self.imap_unordered(fetch_data_worker, [(d, particle_equil_events, self) for d in output_dirs]): #Here we process the returned data from a single directory for state, data in result.items(): if state not in state_data: @@ -555,7 +555,7 @@ def fetch_data(self, particle_equil_events, only_current_statevars = False): #If there's no data, then there's no columns so the next #bit fails. Avoid that if len(df) == 0: - return df + return df, state_data #Here, we're just adjusting the column order to follow what was given by the user. cols = list(df.columns.values) @@ -564,8 +564,6 @@ def fetch_data(self, particle_equil_events, only_current_statevars = False): df = df[[statevar for statevar in self.used_statevariables]+cols] #Now we sort items by the state variables in the order given. df = df.sort_values(by=[statevar for statevar in self.used_statevariables]) + pickle.dump(df, open(self.workdir+".df.pkl", 'wb')) - ##Now we write out the data - - return df - \ No newline at end of file + return df, state_data \ No newline at end of file From 9159dee21f09e30bea5d8be195f12989e8801261 Mon Sep 17 00:00:00 2001 From: Marcus Bannerman Date: Tue, 29 Apr 2025 15:45:05 +0100 Subject: [PATCH 03/19] Updated chunglu plotter --- scripts/plotter.ipynb | 963 +----------------------------------------- 1 file changed, 22 insertions(+), 941 deletions(-) diff --git a/scripts/plotter.ipynb b/scripts/plotter.ipynb index 3b1a19f68..8ce58b075 100644 --- a/scripts/plotter.ipynb +++ b/scripts/plotter.ipynb @@ -2,873 +2,10 @@ "cells": [ { "cell_type": "code", - "execution_count": 53, + "execution_count": null, "id": "bbb1036b", "metadata": {}, - "outputs": [ - { - "data": { - "text/html": [ - "
\n", - "\n", - "\n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - "
NEventsTottTotalpcvu
InitStateLambdakTndensityN
FCC1.51.00.0140004000000023845.0359320.009358+/-0.0000090.148+/-0.004-0.13712+/-0.00019
0.104000400000001247.1170080.0319+/-0.000517+/-6-2.57+/-0.06
0.50400040000000682.734445-0.039+/-0.0043.2+/-0.6-4.648+/-0.013
1.00400040000000352.0543658.353+/-0.0070.284+/-0.010-6.9236+/-0.0008
1.30400040000000153.40526040.61+/-0.050.03097+/-0.00030-8.96942+/-0.00004
1.50.0140004000000025590.7012910.014634+/-0.0000110.0441+/-0.0013-0.09652+/-0.00014
0.104000400000002671.5543740.11652+/-0.000160.426+/-0.020-0.9247+/-0.0009
0.50400040000000655.9841670.5409+/-0.00210.407+/-0.018-3.7988+/-0.0009
1.00400040000000299.92020113.437+/-0.0100.1274+/-0.0020-6.82028+/-0.00023
1.30400040000000121.66455963.43+/-0.060.02037+/-0.00027-8.95679+/-0.00005
2.00.0140004000000025019.7893740.019776+/-0.0000160.0202+/-0.0004-0.08169+/-0.00013
0.104000400000002599.1760140.18250+/-0.000150.171+/-0.005-0.7881+/-0.0005
0.50400040000000574.3013801.3371+/-0.00330.185+/-0.004-3.6622+/-0.0006
1.00400040000000265.20732018.582+/-0.0190.0699+/-0.0020-6.77118+/-0.00028
1.30400040000000103.53887986.59+/-0.100.01361+/-0.00016-8.94862+/-0.00004
2.50.0140004000000023962.1589450.024930+/-0.0000230.01144+/-0.00015-0.07403+/-0.00011
0.104000400000002460.5163920.24633+/-0.000320.0953+/-0.0018-0.7248+/-0.0005
0.50400040000000516.0677962.1419+/-0.00240.114+/-0.004-3.58867+/-0.00035
1.00400040000000240.08835823.714+/-0.0190.0465+/-0.0015-6.74276+/-0.00026
1.3040004000000091.497844110.05+/-0.100.00965+/-0.00009-8.94287+/-0.00006
3.00.0140004000000022859.6378920.030053+/-0.0000230.00768+/-0.00009-0.06928+/-0.00005
0.104000400000002326.2656720.30937+/-0.000280.0581+/-0.0012-0.68636+/-0.00024
0.50400040000000471.6295082.9577+/-0.00320.0778+/-0.0020-3.54119+/-0.00031
1.00400040000000220.99661328.857+/-0.0270.0330+/-0.0007-6.72278+/-0.00026
1.3040004000000082.885040133.37+/-0.130.00727+/-0.00010-8.93878+/-0.00009
inf0.0140004000000048493.5669770.010216+/-0.0000100.0+/-00.0+/-0
0.104000400000004638.7246700.12396+/-0.000110.0+/-00.0+/-0
0.50400040000000814.1862161.6356+/-0.00100.0+/-00.0+/-0
1.00400040000000398.44075510.252+/-0.0100.0+/-00.0+/-0
1.30400040000000136.30320747.68+/-0.040.0+/-00.0+/-0
2.01.00.01400040000000461.9272110.00368+/-0.00021(4.1+/-1.0)e+02-6.0+/-0.4
0.10400040000000198.664502-0.0025+/-0.0010(2.0+/-0.7)e+02-13.17+/-0.29
0.50400040000000143.602432-0.646+/-0.01886+/-20-16.17+/-0.18
1.00400040000000164.077613-28.683+/-0.0072.08+/-0.15-19.9221+/-0.0016
1.30400040000000150.88231547.854+/-0.0340.001190+/-0.000019-21.001139+/-0.000010
1.50.0140004000000015156.3574710.013198+/-0.0000100.1713+/-0.0029-0.2990+/-0.0004
0.10400040000000256.3329240.0033+/-0.000435+/-15-9.76+/-0.16
0.50400040000000173.414606-0.400+/-0.02740+/-17-13.36+/-0.14
1.00400040000000125.644637-25.755+/-0.0191.65+/-0.09-18.8927+/-0.0017
1.30400040000000124.82870571.75+/-0.070.000360+/-0.000004-21.0008050+/-0.0000035
2.00.0140004000000015303.8735260.018525+/-0.0000210.0692+/-0.0016-0.24377+/-0.00024
0.10400040000000448.0477390.0400+/-0.001331+/-10-6.12+/-0.19
0.50400040000000269.417239-0.231+/-0.0178.9+/-2.9-9.67+/-0.08
1.00400040000000111.497971-18.400+/-0.0310.74+/-0.04-18.3158+/-0.0018
1.30400040000000108.68044195.57+/-0.080.0001729+/-0.0000011-21.0006807+/-0.0000027
2.50.0140004000000014762.0702270.023735+/-0.0000130.0371+/-0.0005-0.21868+/-0.00021
0.104000400000001410.9479040.13783+/-0.000300.74+/-0.05-2.204+/-0.004
0.50400040000000311.9337870.238+/-0.0040.263+/-0.019-8.3044+/-0.0013
1.00400040000000101.762522-11.65+/-0.040.448+/-0.015-18.0130+/-0.0015
1.3040004000000097.506480119.37+/-0.100.0001011+/-0.0000007-21.0006130+/-0.0000022
3.00.0140004000000014116.9469080.028891+/-0.0000250.02296+/-0.00028-0.20397+/-0.00016
0.104000400000001426.6386580.2057+/-0.00050.280+/-0.013-1.9688+/-0.0020
0.50400040000000290.4862491.0408+/-0.00320.162+/-0.010-8.2026+/-0.0006
1.0040004000000093.828259-5.50+/-0.040.310+/-0.007-17.8273+/-0.0012
1.3040004000000089.105293143.35+/-0.06(6.44+/-0.04)e-05-21.0005735+/-0.0000031
inf0.0140004000000029776.4998140.0101969+/-0.00000340.0+/-00.0+/-0
0.104000400000002944.7749580.12377+/-0.000120.0+/-00.0+/-0
0.50400040000000532.2841541.6372+/-0.00120.0+/-00.0+/-0
1.00400040000000165.81394010.247+/-0.0110.0+/-00.0+/-0
1.30400040000000155.36996947.78+/-0.040.0+/-00.0+/-0
\n", - "
" - ], - "text/plain": [ - " NEventsTot tTotal \\\n", - "InitState Lambda kT ndensity N \n", - "FCC 1.5 1.0 0.01 4000 40000000 23845.035932 \n", - " 0.10 4000 40000000 1247.117008 \n", - " 0.50 4000 40000000 682.734445 \n", - " 1.00 4000 40000000 352.054365 \n", - " 1.30 4000 40000000 153.405260 \n", - " 1.5 0.01 4000 40000000 25590.701291 \n", - " 0.10 4000 40000000 2671.554374 \n", - " 0.50 4000 40000000 655.984167 \n", - " 1.00 4000 40000000 299.920201 \n", - " 1.30 4000 40000000 121.664559 \n", - " 2.0 0.01 4000 40000000 25019.789374 \n", - " 0.10 4000 40000000 2599.176014 \n", - " 0.50 4000 40000000 574.301380 \n", - " 1.00 4000 40000000 265.207320 \n", - " 1.30 4000 40000000 103.538879 \n", - " 2.5 0.01 4000 40000000 23962.158945 \n", - " 0.10 4000 40000000 2460.516392 \n", - " 0.50 4000 40000000 516.067796 \n", - " 1.00 4000 40000000 240.088358 \n", - " 1.30 4000 40000000 91.497844 \n", - " 3.0 0.01 4000 40000000 22859.637892 \n", - " 0.10 4000 40000000 2326.265672 \n", - " 0.50 4000 40000000 471.629508 \n", - " 1.00 4000 40000000 220.996613 \n", - " 1.30 4000 40000000 82.885040 \n", - " inf 0.01 4000 40000000 48493.566977 \n", - " 0.10 4000 40000000 4638.724670 \n", - " 0.50 4000 40000000 814.186216 \n", - " 1.00 4000 40000000 398.440755 \n", - " 1.30 4000 40000000 136.303207 \n", - " 2.0 1.0 0.01 4000 40000000 461.927211 \n", - " 0.10 4000 40000000 198.664502 \n", - " 0.50 4000 40000000 143.602432 \n", - " 1.00 4000 40000000 164.077613 \n", - " 1.30 4000 40000000 150.882315 \n", - " 1.5 0.01 4000 40000000 15156.357471 \n", - " 0.10 4000 40000000 256.332924 \n", - " 0.50 4000 40000000 173.414606 \n", - " 1.00 4000 40000000 125.644637 \n", - " 1.30 4000 40000000 124.828705 \n", - " 2.0 0.01 4000 40000000 15303.873526 \n", - " 0.10 4000 40000000 448.047739 \n", - " 0.50 4000 40000000 269.417239 \n", - " 1.00 4000 40000000 111.497971 \n", - " 1.30 4000 40000000 108.680441 \n", - " 2.5 0.01 4000 40000000 14762.070227 \n", - " 0.10 4000 40000000 1410.947904 \n", - " 0.50 4000 40000000 311.933787 \n", - " 1.00 4000 40000000 101.762522 \n", - " 1.30 4000 40000000 97.506480 \n", - " 3.0 0.01 4000 40000000 14116.946908 \n", - " 0.10 4000 40000000 1426.638658 \n", - " 0.50 4000 40000000 290.486249 \n", - " 1.00 4000 40000000 93.828259 \n", - " 1.30 4000 40000000 89.105293 \n", - " inf 0.01 4000 40000000 29776.499814 \n", - " 0.10 4000 40000000 2944.774958 \n", - " 0.50 4000 40000000 532.284154 \n", - " 1.00 4000 40000000 165.813940 \n", - " 1.30 4000 40000000 155.369969 \n", - "\n", - " p \\\n", - "InitState Lambda kT ndensity N \n", - "FCC 1.5 1.0 0.01 4000 0.009358+/-0.000009 \n", - " 0.10 4000 0.0319+/-0.0005 \n", - " 0.50 4000 -0.039+/-0.004 \n", - " 1.00 4000 8.353+/-0.007 \n", - " 1.30 4000 40.61+/-0.05 \n", - " 1.5 0.01 4000 0.014634+/-0.000011 \n", - " 0.10 4000 0.11652+/-0.00016 \n", - " 0.50 4000 0.5409+/-0.0021 \n", - " 1.00 4000 13.437+/-0.010 \n", - " 1.30 4000 63.43+/-0.06 \n", - " 2.0 0.01 4000 0.019776+/-0.000016 \n", - " 0.10 4000 0.18250+/-0.00015 \n", - " 0.50 4000 1.3371+/-0.0033 \n", - " 1.00 4000 18.582+/-0.019 \n", - " 1.30 4000 86.59+/-0.10 \n", - " 2.5 0.01 4000 0.024930+/-0.000023 \n", - " 0.10 4000 0.24633+/-0.00032 \n", - " 0.50 4000 2.1419+/-0.0024 \n", - " 1.00 4000 23.714+/-0.019 \n", - " 1.30 4000 110.05+/-0.10 \n", - " 3.0 0.01 4000 0.030053+/-0.000023 \n", - " 0.10 4000 0.30937+/-0.00028 \n", - " 0.50 4000 2.9577+/-0.0032 \n", - " 1.00 4000 28.857+/-0.027 \n", - " 1.30 4000 133.37+/-0.13 \n", - " inf 0.01 4000 0.010216+/-0.000010 \n", - " 0.10 4000 0.12396+/-0.00011 \n", - " 0.50 4000 1.6356+/-0.0010 \n", - " 1.00 4000 10.252+/-0.010 \n", - " 1.30 4000 47.68+/-0.04 \n", - " 2.0 1.0 0.01 4000 0.00368+/-0.00021 \n", - " 0.10 4000 -0.0025+/-0.0010 \n", - " 0.50 4000 -0.646+/-0.018 \n", - " 1.00 4000 -28.683+/-0.007 \n", - " 1.30 4000 47.854+/-0.034 \n", - " 1.5 0.01 4000 0.013198+/-0.000010 \n", - " 0.10 4000 0.0033+/-0.0004 \n", - " 0.50 4000 -0.400+/-0.027 \n", - " 1.00 4000 -25.755+/-0.019 \n", - " 1.30 4000 71.75+/-0.07 \n", - " 2.0 0.01 4000 0.018525+/-0.000021 \n", - " 0.10 4000 0.0400+/-0.0013 \n", - " 0.50 4000 -0.231+/-0.017 \n", - " 1.00 4000 -18.400+/-0.031 \n", - " 1.30 4000 95.57+/-0.08 \n", - " 2.5 0.01 4000 0.023735+/-0.000013 \n", - " 0.10 4000 0.13783+/-0.00030 \n", - " 0.50 4000 0.238+/-0.004 \n", - " 1.00 4000 -11.65+/-0.04 \n", - " 1.30 4000 119.37+/-0.10 \n", - " 3.0 0.01 4000 0.028891+/-0.000025 \n", - " 0.10 4000 0.2057+/-0.0005 \n", - " 0.50 4000 1.0408+/-0.0032 \n", - " 1.00 4000 -5.50+/-0.04 \n", - " 1.30 4000 143.35+/-0.06 \n", - " inf 0.01 4000 0.0101969+/-0.0000034 \n", - " 0.10 4000 0.12377+/-0.00012 \n", - " 0.50 4000 1.6372+/-0.0012 \n", - " 1.00 4000 10.247+/-0.011 \n", - " 1.30 4000 47.78+/-0.04 \n", - "\n", - " cv \\\n", - "InitState Lambda kT ndensity N \n", - "FCC 1.5 1.0 0.01 4000 0.148+/-0.004 \n", - " 0.10 4000 17+/-6 \n", - " 0.50 4000 3.2+/-0.6 \n", - " 1.00 4000 0.284+/-0.010 \n", - " 1.30 4000 0.03097+/-0.00030 \n", - " 1.5 0.01 4000 0.0441+/-0.0013 \n", - " 0.10 4000 0.426+/-0.020 \n", - " 0.50 4000 0.407+/-0.018 \n", - " 1.00 4000 0.1274+/-0.0020 \n", - " 1.30 4000 0.02037+/-0.00027 \n", - " 2.0 0.01 4000 0.0202+/-0.0004 \n", - " 0.10 4000 0.171+/-0.005 \n", - " 0.50 4000 0.185+/-0.004 \n", - " 1.00 4000 0.0699+/-0.0020 \n", - " 1.30 4000 0.01361+/-0.00016 \n", - " 2.5 0.01 4000 0.01144+/-0.00015 \n", - " 0.10 4000 0.0953+/-0.0018 \n", - " 0.50 4000 0.114+/-0.004 \n", - " 1.00 4000 0.0465+/-0.0015 \n", - " 1.30 4000 0.00965+/-0.00009 \n", - " 3.0 0.01 4000 0.00768+/-0.00009 \n", - " 0.10 4000 0.0581+/-0.0012 \n", - " 0.50 4000 0.0778+/-0.0020 \n", - " 1.00 4000 0.0330+/-0.0007 \n", - " 1.30 4000 0.00727+/-0.00010 \n", - " inf 0.01 4000 0.0+/-0 \n", - " 0.10 4000 0.0+/-0 \n", - " 0.50 4000 0.0+/-0 \n", - " 1.00 4000 0.0+/-0 \n", - " 1.30 4000 0.0+/-0 \n", - " 2.0 1.0 0.01 4000 (4.1+/-1.0)e+02 \n", - " 0.10 4000 (2.0+/-0.7)e+02 \n", - " 0.50 4000 86+/-20 \n", - " 1.00 4000 2.08+/-0.15 \n", - " 1.30 4000 0.001190+/-0.000019 \n", - " 1.5 0.01 4000 0.1713+/-0.0029 \n", - " 0.10 4000 35+/-15 \n", - " 0.50 4000 40+/-17 \n", - " 1.00 4000 1.65+/-0.09 \n", - " 1.30 4000 0.000360+/-0.000004 \n", - " 2.0 0.01 4000 0.0692+/-0.0016 \n", - " 0.10 4000 31+/-10 \n", - " 0.50 4000 8.9+/-2.9 \n", - " 1.00 4000 0.74+/-0.04 \n", - " 1.30 4000 0.0001729+/-0.0000011 \n", - " 2.5 0.01 4000 0.0371+/-0.0005 \n", - " 0.10 4000 0.74+/-0.05 \n", - " 0.50 4000 0.263+/-0.019 \n", - " 1.00 4000 0.448+/-0.015 \n", - " 1.30 4000 0.0001011+/-0.0000007 \n", - " 3.0 0.01 4000 0.02296+/-0.00028 \n", - " 0.10 4000 0.280+/-0.013 \n", - " 0.50 4000 0.162+/-0.010 \n", - " 1.00 4000 0.310+/-0.007 \n", - " 1.30 4000 (6.44+/-0.04)e-05 \n", - " inf 0.01 4000 0.0+/-0 \n", - " 0.10 4000 0.0+/-0 \n", - " 0.50 4000 0.0+/-0 \n", - " 1.00 4000 0.0+/-0 \n", - " 1.30 4000 0.0+/-0 \n", - "\n", - " u \n", - "InitState Lambda kT ndensity N \n", - "FCC 1.5 1.0 0.01 4000 -0.13712+/-0.00019 \n", - " 0.10 4000 -2.57+/-0.06 \n", - " 0.50 4000 -4.648+/-0.013 \n", - " 1.00 4000 -6.9236+/-0.0008 \n", - " 1.30 4000 -8.96942+/-0.00004 \n", - " 1.5 0.01 4000 -0.09652+/-0.00014 \n", - " 0.10 4000 -0.9247+/-0.0009 \n", - " 0.50 4000 -3.7988+/-0.0009 \n", - " 1.00 4000 -6.82028+/-0.00023 \n", - " 1.30 4000 -8.95679+/-0.00005 \n", - " 2.0 0.01 4000 -0.08169+/-0.00013 \n", - " 0.10 4000 -0.7881+/-0.0005 \n", - " 0.50 4000 -3.6622+/-0.0006 \n", - " 1.00 4000 -6.77118+/-0.00028 \n", - " 1.30 4000 -8.94862+/-0.00004 \n", - " 2.5 0.01 4000 -0.07403+/-0.00011 \n", - " 0.10 4000 -0.7248+/-0.0005 \n", - " 0.50 4000 -3.58867+/-0.00035 \n", - " 1.00 4000 -6.74276+/-0.00026 \n", - " 1.30 4000 -8.94287+/-0.00006 \n", - " 3.0 0.01 4000 -0.06928+/-0.00005 \n", - " 0.10 4000 -0.68636+/-0.00024 \n", - " 0.50 4000 -3.54119+/-0.00031 \n", - " 1.00 4000 -6.72278+/-0.00026 \n", - " 1.30 4000 -8.93878+/-0.00009 \n", - " inf 0.01 4000 0.0+/-0 \n", - " 0.10 4000 0.0+/-0 \n", - " 0.50 4000 0.0+/-0 \n", - " 1.00 4000 0.0+/-0 \n", - " 1.30 4000 0.0+/-0 \n", - " 2.0 1.0 0.01 4000 -6.0+/-0.4 \n", - " 0.10 4000 -13.17+/-0.29 \n", - " 0.50 4000 -16.17+/-0.18 \n", - " 1.00 4000 -19.9221+/-0.0016 \n", - " 1.30 4000 -21.001139+/-0.000010 \n", - " 1.5 0.01 4000 -0.2990+/-0.0004 \n", - " 0.10 4000 -9.76+/-0.16 \n", - " 0.50 4000 -13.36+/-0.14 \n", - " 1.00 4000 -18.8927+/-0.0017 \n", - " 1.30 4000 -21.0008050+/-0.0000035 \n", - " 2.0 0.01 4000 -0.24377+/-0.00024 \n", - " 0.10 4000 -6.12+/-0.19 \n", - " 0.50 4000 -9.67+/-0.08 \n", - " 1.00 4000 -18.3158+/-0.0018 \n", - " 1.30 4000 -21.0006807+/-0.0000027 \n", - " 2.5 0.01 4000 -0.21868+/-0.00021 \n", - " 0.10 4000 -2.204+/-0.004 \n", - " 0.50 4000 -8.3044+/-0.0013 \n", - " 1.00 4000 -18.0130+/-0.0015 \n", - " 1.30 4000 -21.0006130+/-0.0000022 \n", - " 3.0 0.01 4000 -0.20397+/-0.00016 \n", - " 0.10 4000 -1.9688+/-0.0020 \n", - " 0.50 4000 -8.2026+/-0.0006 \n", - " 1.00 4000 -17.8273+/-0.0012 \n", - " 1.30 4000 -21.0005735+/-0.0000031 \n", - " inf 0.01 4000 0.0+/-0 \n", - " 0.10 4000 0.0+/-0 \n", - " 0.50 4000 0.0+/-0 \n", - " 1.00 4000 0.0+/-0 \n", - " 1.30 4000 0.0+/-0 " - ] - }, - "execution_count": 53, - "metadata": {}, - "output_type": "execute_result" - } - ], + "outputs": [], "source": [ "import pandas\n", "import pickle\n", @@ -894,7 +31,6 @@ "# Filter out columns that do not change (boring state variables i.e.)\n", "#orig_df = orig_df.loc[:, orig_df.nunique() > 1]\n", "\n", - "\n", "# Create a dataframe of scalars\n", "scalar_df = orig_df.copy()\n", "for col in scalar_df.columns:\n", @@ -906,84 +42,29 @@ }, { "cell_type": "code", - "execution_count": 55, + "execution_count": null, "id": "3c5f65e0", "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "InitState Lambda kT ndensity N \n", - "FCC 1.5 1.0 0.01 4000 [(1, 1), (2, 2), (1, 2), (1, 3), (4, 5), (3, 5...\n", - " 0.10 4000 [(8, 10), (10, 10), (6, 10), (10, 11), (7, 10)...\n", - " 0.50 4000 [(8, 9), (8, 10), (8, 8), (8, 11), (8, 12), (4...\n", - " 1.00 4000 [(13, 16), (13, 14), (13, 13), (12, 13), (14, ...\n", - " 1.30 4000 [(18, 18), (17, 18), (17, 17), (16, 18), (16, ...\n", - " 1.5 0.01 4000 [(1, 1), (2, 3), (1, 3), (1, 2), (2, 2), (2, 4...\n", - " 0.10 4000 [(1, 2), (1, 1), (2, 4), (2, 3), (3, 3), (3, 5...\n", - " 0.50 4000 [(6, 7), (6, 8), (5, 6), (6, 9), (7, 8), (8, 1...\n", - " 1.00 4000 [(12, 14), (13, 14), (14, 14), (14, 15), (13, ...\n", - " 1.30 4000 [(18, 18), (17, 18), (16, 18), (17, 17), (16, ...\n", - " 2.0 0.01 4000 [(1, 1), (2, 2), (1, 2), (1, 3), (2, 3), (3, 3)]\n", - " 0.10 4000 [(1, 3), (2, 4), (2, 3), (2, 2), (4, 6), (4, 4...\n", - " 0.50 4000 [(7, 8), (8, 10), (8, 9), (8, 8), (6, 8), (6, ...\n", - " 1.00 4000 [(14, 14), (14, 15), (13, 14), (12, 14), (13, ...\n", - " 1.30 4000 [(17, 18), (18, 18), (16, 18), (17, 17), (16, ...\n", - " 2.5 0.01 4000 [(1, 1), (1, 2), (2, 2), (1, 3), (2, 3), (3, 3)]\n", - " 0.10 4000 [(1, 2), (1, 3), (2, 3), (1, 1), (3, 4), (4, 4...\n", - " 0.50 4000 [(5, 5), (5, 7), (5, 8), (8, 8), (8, 10), (7, ...\n", - " 1.00 4000 [(14, 14), (13, 14), (12, 14), (14, 15), (12, ...\n", - " 1.30 4000 [(18, 18), (17, 18), (16, 18), (17, 17), (16, ...\n", - " 3.0 0.01 4000 [(1, 2), (1, 1), (2, 2), (1, 3), (2, 3), (3, 3)]\n", - " 0.10 4000 [(2, 3), (1, 1), (1, 2), (2, 4), (4, 4), (1, 3...\n", - " 0.50 4000 [(9, 10), (8, 9), (9, 11), (7, 9), (6, 9), (9,...\n", - " 1.00 4000 [(12, 13), (12, 14), (12, 15), (13, 14), (14, ...\n", - " 1.30 4000 [(18, 18), (17, 18), (17, 17), (16, 18), (16, ...\n", - " inf 0.01 4000 [(1, 1), (1, 2), (2, 2), (2, 3), (1, 3), (3, 3)]\n", - " 0.10 4000 [(2, 2), (1, 2), (2, 3), (1, 3), (1, 4), (1, 1...\n", - " 0.50 4000 [(6, 6), (5, 6), (6, 10), (6, 8), (5, 7), (6, ...\n", - " 1.00 4000 [(13, 15), (13, 13), (13, 14), (11, 13), (12, ...\n", - " 1.30 4000 [(18, 18), (17, 18), (16, 18), (17, 17), (16, ...\n", - " 2.0 1.0 0.01 4000 [(21, 32), (21, 48), (21, 41), (21, 34), (21, ...\n", - " 0.10 4000 [(38, 38), (37, 38), (22, 38), (24, 38), (19, ...\n", - " 0.50 4000 [(34, 45), (45, 45), (43, 45), (36, 45), (22, ...\n", - " 1.00 4000 [(39, 40), (40, 41), (40, 42), (40, 40), (38, ...\n", - " 1.30 4000 [(42, 42), (42, 43), (43, 43)]\n", - " 1.5 0.01 4000 [(1, 1), (1, 2), (2, 2), (2, 3), (4, 4), (3, 3...\n", - " 0.10 4000 [(12, 14), (12, 13), (12, 20), (12, 19), (12, ...\n", - " 0.50 4000 [(29, 31), (31, 32), (31, 31), (24, 31), (27, ...\n", - " 1.00 4000 [(38, 38), (33, 38), (36, 38), (35, 38), (38, ...\n", - " 1.30 4000 [(42, 42), (42, 43), (43, 43)]\n", - " 2.0 0.01 4000 [(1, 1), (2, 3), (1, 3), (3, 3), (2, 2), (1, 2...\n", - " 0.10 4000 [(2, 2), (22, 22), (21, 22), (22, 24), (15, 22...\n", - " 0.50 4000 [(18, 21), (17, 18), (18, 19), (18, 20), (16, ...\n", - " 1.00 4000 [(34, 41), (34, 39), (34, 37), (34, 38), (34, ...\n", - " 1.30 4000 [(42, 42), (42, 43), (43, 43)]\n", - " 2.5 0.01 4000 [(1, 2), (1, 3), (1, 1), (2, 2), (2, 4), (3, 4...\n", - " 0.10 4000 [(4, 8), (4, 6), (3, 4), (6, 11), (5, 6), (6, ...\n", - " 0.50 4000 [(18, 20), (17, 20), (20, 21), (19, 20), (20, ...\n", - " 1.00 4000 [(31, 37), (37, 38), (35, 37), (36, 37), (34, ...\n", - " 1.30 4000 [(42, 42), (42, 43), (43, 43)]\n", - " 3.0 0.01 4000 [(1, 3), (1, 1), (1, 2), (2, 3), (2, 2), (3, 3...\n", - " 0.10 4000 [(5, 8), (8, 8), (6, 8), (7, 8), (1, 3), (2, 6...\n", - " 0.50 4000 [(16, 16), (15, 16), (16, 17), (16, 18), (16, ...\n", - " 1.00 4000 [(36, 38), (35, 38), (38, 39), (38, 38), (38, ...\n", - " 1.30 4000 [(42, 42), (42, 43), (43, 43)]\n", - " inf 0.01 4000 [(1, 1), (1, 3), (2, 2), (2, 3), (1, 2), (3, 3...\n", - " 0.10 4000 [(3, 4), (3, 3), (2, 3), (5, 6), (6, 7), (3, 6...\n", - " 0.50 4000 [(13, 16), (15, 16), (16, 17), (14, 16), (16, ...\n", - " 1.00 4000 [(36, 37), (36, 36), (35, 36), (31, 36), (30, ...\n", - " 1.30 4000 [(42, 42), (42, 43), (43, 43)]\n", - "Name: ChungLu, dtype: object" - ] - }, - "execution_count": 55, - "metadata": {}, - "output_type": "execute_result" - } - ], + "outputs": [], "source": [ - "orig_df[\"ChungLu\"]" + "import matplotlib.pyplot as plt\n", + "\n", + "for state in state_data.keys():\n", + " if not (\"Lambda\", 1.5) in state or not (\"kT\", 1.5) in state or not (\"ndensity\", 1.0) in state:\n", + " continue\n", + "\n", + " row = state_data[state][\"ChungLu\"]\n", + " curve = row.ufloat()\n", + " curve = [(k[0] * k[1], v.nominal_value, v.std_dev) for k,v in curve.items()]\n", + " print(curve)\n", + " plt.errorbar(\n", + " [k[0] for k in curve],\n", + " [k[1] for k in curve],\n", + " yerr=[k[2] for k in curve],\n", + " label=state[-1]\n", + " )\n", + "plt.legend()\n", + "plt.show() " ] }, { From 7939c107a15b1fed564c1f8e41b879f974ee63e6 Mon Sep 17 00:00:00 2001 From: Marcus Bannerman Date: Wed, 30 Apr 2025 09:51:00 +0100 Subject: [PATCH 04/19] Removing windows builds for now --- src/pydynamo/test_weighted_types.py | 8 ++++---- src/pydynamo/weighted_types.py | 8 +++----- 2 files changed, 7 insertions(+), 9 deletions(-) diff --git a/src/pydynamo/test_weighted_types.py b/src/pydynamo/test_weighted_types.py index 6f6f25335..debd7670b 100644 --- a/src/pydynamo/test_weighted_types.py +++ b/src/pydynamo/test_weighted_types.py @@ -73,10 +73,10 @@ def test_keyed_array(): assert c["b"] == pytest.approx(3) assert c["c"] == pytest.approx(4.5) - c = a / b - assert c["a"] == pytest.approx(0.5) - assert c["b"] == pytest.approx(0.5) - assert "c" not in c + # Check division of two keyed arrays fails + with pytest.raises(Exception): + c = a / b + c = a / 2 assert c["a"] == pytest.approx(0.5) assert c["b"] == pytest.approx(1) diff --git a/src/pydynamo/weighted_types.py b/src/pydynamo/weighted_types.py index 86ec02b92..2a2aca923 100644 --- a/src/pydynamo/weighted_types.py +++ b/src/pydynamo/weighted_types.py @@ -10,11 +10,9 @@ class KeyedArray(): """A key-value store of values that has element-wise addition and multiplication. Any missing values are assumed to be zero. This is needed for observations that do not appear in some simulations. """ - - store = defaultdict(float) - - def __init__(self, values = None): - self.store = defaultdict(float) + def __init__(self, values = None, type = float): + self.type = type + self.store = defaultdict(type) if values is not None: if isinstance(values, KeyedArray): for k, v in values.items(): From 062f66fd541306df38ce5f9038013494d0427627 Mon Sep 17 00:00:00 2001 From: Marcus Bannerman Date: Wed, 30 Apr 2025 10:22:33 +0100 Subject: [PATCH 05/19] Made python builds part of the CI --- src/pydynamo/test_weighted_types.py | 27 ++++++++++++++++++++++++++- 1 file changed, 26 insertions(+), 1 deletion(-) diff --git a/src/pydynamo/test_weighted_types.py b/src/pydynamo/test_weighted_types.py index debd7670b..a55e62a05 100644 --- a/src/pydynamo/test_weighted_types.py +++ b/src/pydynamo/test_weighted_types.py @@ -96,4 +96,29 @@ def test_keyed_array(): assert c.avg()["b"] == pytest.approx(c2.avg()) for i in range(3): assert c.stats()[i]["a"] == pytest.approx(c1.stats()[i]) - assert c.stats()[i]["b"] == pytest.approx(c2.stats()[i]) \ No newline at end of file + assert c.stats()[i]["b"] == pytest.approx(c2.stats()[i]) + +def test_keyed_keyed_array(): + a = KeyedArray(type=KeyedArray) + a["a"] = KeyedArray({"b":1}) + a["a"] = KeyedArray({"c":3}) + a["b"] = KeyedArray({"c":3}) + a["c"] = KeyedArray({"d":4}) + + b = KeyedArray(type=KeyedArray) + b["a"] = KeyedArray({"b":2}) + + c = a + b + assert c["a"]["b"] == pytest.approx(3) + assert c["b"]["c"] == pytest.approx(3) + assert c["c"]["d"] == pytest.approx(4) + + c = a - b + assert c["a"]["b"] == pytest.approx(-1) + assert c["b"]["c"] == pytest.approx(3) + assert c["c"]["d"] == pytest.approx(4) + + c = a * b + assert c["a"]["b"] == pytest.approx(2) + assert "b" not in c + assert "c" not in c \ No newline at end of file From bb78b01c6559a88016dc11e5c5a4891ba01c1bd4 Mon Sep 17 00:00:00 2001 From: Marcus Bannerman Date: Wed, 30 Apr 2025 15:30:01 +0100 Subject: [PATCH 06/19] Extended for nested arrays --- src/pydynamo/test_weighted_types.py | 33 ++++++++++++++++++----------- src/pydynamo/weighted_types.py | 19 ++++++++--------- 2 files changed, 30 insertions(+), 22 deletions(-) diff --git a/src/pydynamo/test_weighted_types.py b/src/pydynamo/test_weighted_types.py index a55e62a05..b1ea9802d 100644 --- a/src/pydynamo/test_weighted_types.py +++ b/src/pydynamo/test_weighted_types.py @@ -100,25 +100,34 @@ def test_keyed_array(): def test_keyed_keyed_array(): a = KeyedArray(type=KeyedArray) - a["a"] = KeyedArray({"b":1}) - a["a"] = KeyedArray({"c":3}) - a["b"] = KeyedArray({"c":3}) - a["c"] = KeyedArray({"d":4}) + a["a"]["1"] = 1 + a["a"]["2"] = 3 + a["b"]["2"] = 3 + a["c"]["3"] = 4 + + # Check items that are there + assert a["a"]["1"] == 1 + assert a["a"]["2"] == 3 + assert a["b"]["2"] == 3 + assert a["c"]["3"] == 4 + + # Check items that are not there + assert a["a"]["3"] == 0 b = KeyedArray(type=KeyedArray) - b["a"] = KeyedArray({"b":2}) + b["a"] = KeyedArray(float, {"1":2}) c = a + b - assert c["a"]["b"] == pytest.approx(3) - assert c["b"]["c"] == pytest.approx(3) - assert c["c"]["d"] == pytest.approx(4) + assert c["a"]["1"] == pytest.approx(3) + assert c["b"]["2"] == pytest.approx(3) + assert c["c"]["3"] == pytest.approx(4) c = a - b - assert c["a"]["b"] == pytest.approx(-1) - assert c["b"]["c"] == pytest.approx(3) - assert c["c"]["d"] == pytest.approx(4) + assert c["a"]["1"] == pytest.approx(-1) + assert c["b"]["2"] == pytest.approx(3) + assert c["c"]["3"] == pytest.approx(4) c = a * b - assert c["a"]["b"] == pytest.approx(2) + assert c["a"]["1"] == pytest.approx(2) assert "b" not in c assert "c" not in c \ No newline at end of file diff --git a/src/pydynamo/weighted_types.py b/src/pydynamo/weighted_types.py index 2a2aca923..392da5806 100644 --- a/src/pydynamo/weighted_types.py +++ b/src/pydynamo/weighted_types.py @@ -10,16 +10,15 @@ class KeyedArray(): """A key-value store of values that has element-wise addition and multiplication. Any missing values are assumed to be zero. This is needed for observations that do not appear in some simulations. """ - def __init__(self, values = None, type = float): + def __init__(self, type = float, values = None): self.type = type + #print("KeyedArray(", repr(type),",", repr(values),")") self.store = defaultdict(type) + import copy if values is not None: - if isinstance(values, KeyedArray): - for k, v in values.items(): - self.store[k] = v - elif isinstance(values, dict): + if isinstance(values, KeyedArray) or isinstance(values, dict): for k, v in values.items(): - self.store[k] = v + self.store[k] = copy.copy(v) else: raise RuntimeError("Cannot create KeyedArray from non-dict or non-KeyedArray") @@ -114,10 +113,10 @@ def __truediv__(self, rhs): return retval def __repr__(self): - return str(self.store) + return "KeyedArray("+self.type.__name__ +", {" + ", ".join([repr(k) + ": " + repr(v) for k, v in self.store.items()]) + "})" def __str__(self): - return str(self.store) + return self.__repr__() def element_wise_multiply(a, b): if isinstance(a, numpy.ndarray) or isinstance(b, numpy.ndarray): @@ -135,13 +134,13 @@ def zeros_like(a): def maximum(a, b): if isinstance(a, KeyedArray) and isinstance(b, KeyedArray): - return KeyedArray({k: maximum(a[k], b[k]) for k in set(a.keys()).union(b.keys())}) + return KeyedArray(type=a.type, values = {k: maximum(a[k], b[k]) for k in set(a.keys()).union(b.keys())}) else: return numpy.maximum(a,b) def sqrt(a): if isinstance(a, KeyedArray): - return KeyedArray({k: sqrt(a[k]) for k in a.keys()}) + return KeyedArray(type = a.type, values = {k: sqrt(a[k]) for k in a.keys()}) else: return numpy.sqrt(a) From 710434afb0e7fc294777a76bd9b0ffe9e4c4269c Mon Sep 17 00:00:00 2001 From: Marcus Bannerman Date: Thu, 1 May 2025 13:16:10 +0100 Subject: [PATCH 07/19] Got multiple outputs working from the ChungLu plugin as a template for others in pydynamo --- .vscode/launch.json | 8 ++++ scripts/plotter.ipynb | 17 +++++--- src/pydynamo/__init__.py | 2 +- src/pydynamo/output_properties.py | 68 ++++++++++++++++------------- src/pydynamo/test_weighted_types.py | 21 ++++++++- src/pydynamo/weighted_types.py | 20 +++++++-- 6 files changed, 95 insertions(+), 41 deletions(-) diff --git a/.vscode/launch.json b/.vscode/launch.json index 6afb5f936..07cce0a42 100644 --- a/.vscode/launch.json +++ b/.vscode/launch.json @@ -4,6 +4,14 @@ // For more information, visit: https://go.microsoft.com/fwlink/?linkid=830387 "version": "0.2.0", "configurations": [ + { + "name": "Python Debugger: Current File", + "type": "debugpy", + "request": "launch", + "program": "${file}", + "console": "integratedTerminal", + "cwd": "${workspaceFolder}/scripts" + }, { "name": "(ctest) Launch", "type": "cppdbg", diff --git a/scripts/plotter.ipynb b/scripts/plotter.ipynb index 8ce58b075..338ff45d9 100644 --- a/scripts/plotter.ipynb +++ b/scripts/plotter.ipynb @@ -54,23 +54,30 @@ " continue\n", "\n", " row = state_data[state][\"ChungLu\"]\n", - " curve = row.ufloat()\n", - " curve = [(k[0] * k[1], v.nominal_value, v.std_dev) for k,v in curve.items()]\n", + "\n", + " bond_order = row[\"bond_order_count\"].ufloat()\n", + " particle_order = row[\"order_count\"].ufloat()\n", + "\n", + " curve = [(k[0] * k[1], (v / (particle_order[k[0]] * (particle_order[k[1]] - float(k[0] == k[1])))).nominal_value, (v / (particle_order[k[0]] * (particle_order[k[1]] - float(k[0] == k[1])))).std_dev) for k,v in bond_order.items()]\n", " print(curve)\n", " plt.errorbar(\n", " [k[0] for k in curve],\n", " [k[1] for k in curve],\n", " yerr=[k[2] for k in curve],\n", - " label=state[-1]\n", + " label=state[-1],\n", + " fmt=\"x\",\n", " )\n", + "plt.xlabel(\"$k_i\\\\,k_j$\")\n", + "plt.ylabel(\"$\\\\frac{N(k_i,\\\\,k_j)}{N(k_i)\\\\,N(k_j)}$\")\n", "plt.legend()\n", - "plt.show() " + "plt.grid()\n", + "plt.show()" ] }, { "cell_type": "code", "execution_count": null, - "id": "acb5456f", + "id": "fc6e9d0c", "metadata": {}, "outputs": [], "source": [] diff --git a/src/pydynamo/__init__.py b/src/pydynamo/__init__.py index e3fa0fc25..13a261957 100755 --- a/src/pydynamo/__init__.py +++ b/src/pydynamo/__init__.py @@ -313,7 +313,7 @@ def getnextstatedir(self, state, oldpath = None): return newpath idx += 1 - def imap_unordered(self, func, iterable, chunksize=2): + def imap_unordered(self, func, iterable, chunksize=1): """A version of imap_unordered that works with multiprocessing.Pool""" # Chunksize halves overhead for tiny tasks, but also doesn't limit # parallelism when doing a few slow tasks, or on a system with many diff --git a/src/pydynamo/output_properties.py b/src/pydynamo/output_properties.py index 7d10d2f63..dec71db1d 100644 --- a/src/pydynamo/output_properties.py +++ b/src/pydynamo/output_properties.py @@ -7,7 +7,8 @@ from pydynamo.config_files import ConfigFile from pydynamo.file_types import XMLFile, validate_xmlfile -from pydynamo.weighted_types import KeyedArray, WeightedType +from pydynamo.weighted_types import (KeyedArray, KeyedWeightedKeyedArray, + WeightedType) # A XMLFile/ElementTree but specialised for DynamO output files @@ -106,6 +107,22 @@ def weight(self, outputfile): def result(self, state, outputfile, configfilename, counter, manager, output_dir): return WeightedType(self.value(outputfile), self.weight(outputfile)) +OutputFile.output_props["N"] = SingleAttrib('ParticleCount', 'val', [], [], [], missing_val=None)#We use missing_val=None to cause an error if the tag is missing +OutputFile.output_props["p"] = SingleAttrib('Pressure', 'Avg', [], [], [], missing_val=None) +OutputFile.output_props["cv"] = SingleAttrib('ResidualHeatCapacity', 'Value', [], [], [], div_by_N=True, missing_val=None) +OutputFile.output_props["u"] = SingleAttrib('UConfigurational', 'Mean', [], [], [], div_by_N=True, missing_val=None) +OutputFile.output_props["T"] = SingleAttrib('Temperature', 'Mean', [], [], [], missing_val=None) +OutputFile.output_props["density"] = SingleAttrib('Density', 'val', [], [], [], missing_val=None) +OutputFile.output_props["MSD"] = SingleAttrib('MSD/Species', 'diffusionCoeff', [], [], ['-LMSD'], missing_val=None, skip_missing=True) +OutputFile.output_props["NeventsSO"] = SingleAttrib('EventCounters/Entry[@Name="SOCells"]', # Outputfile tag name + 'Count', # Outputfile tag attribute name + ["Rso"], # Required state variable + [], # Required output variables + [], # Required output plugins + div_by_N=True, # Divide the count by N + div_by_t=True, # Also divide by t + missing_val=0) # If counter is missing, return 0 + def parseToArray(text): data = [] for row in text.split('\n'): @@ -120,6 +137,7 @@ def __init__(self): def result(self, state, outputfile, configfilename, counter, manager, output_dir): return None +OutputFile.output_props["CollisionMatrix"] = CollisionMatrixOutputProperty() class VACFOutputProperty(OutputProperty): def __init__(self): @@ -135,6 +153,7 @@ def result(self, state, outputfile, configfilename, counter, manager, output_dir for tag in outputfile.tree.findall('.//VACF/Topology/Structure'): pickle.dump(parseToArray(tag.text), open(filename_root + '/topology_'+tag.attrib['Name']+'.pkl', 'wb')) return None +OutputFile.output_props["VACF"] = VACFOutputProperty() class RadialDistributionOutputProperty(OutputProperty): def __init__(self): @@ -182,6 +201,7 @@ def result(self, state, outputfile, configfilename, counter, manager, output_dir central_moments[n-1] = sum(scipy.special.comb(n, i) * moments[i] * (N0-central_moments[0,:])**(n-i) for i in range(n+1)) return WeightedType(central_moments, samples) +OutputFile.output_props["RadialDistribution"] = RadialDistributionOutputProperty() class RadialDistEndOutputProperty(OutputProperty): def __init__(self): @@ -209,6 +229,7 @@ def result(self, state, outputfile, configfilename, counter, manager, output_dir output_pkl = filename_root + '/species_'+A+'_'+B+'.pkl' pickle.dump(parseToArray(tag.text), open(output_pkl, 'wb')) return None +OutputFile.output_props["RadialDistEnd"] = RadialDistEndOutputProperty() class OrderParameterProperty(OutputProperty): ''' @@ -250,6 +271,7 @@ def result(self, state, outputfile, configfilename, counter, manager, output_dir ql_value = ql.particle_order return WeightedType(numpy.mean(ql_value), 1) +OutputFile.output_props["FCCOrder"] = OrderParameterProperty(6) class ChungLuConfigurationModel(OutputProperty): ''' @@ -261,7 +283,7 @@ class ChungLuConfigurationModel(OutputProperty): def __init__(self): OutputProperty.__init__(self, dependent_statevars=[], dependent_outputs=[], dependent_outputplugins=[]) - def result(self, state, outputfile, configfilename, counter, manager, output_dir): + def result(self, state, outputfile, configfilename, bond_order_counter, manager, output_dir): configfile = ConfigFile(configfilename) if len(configfile.tree.findall('.//Interaction/CaptureMap')) != 1: @@ -275,47 +297,33 @@ def result(self, state, outputfile, configfilename, counter, manager, output_dir for pair in configfile.tree.findall('.//Interaction/CaptureMap/Pair'): G.add_edge(int(pair.attrib['ID1']), int(pair.attrib['ID2'])) - degrees = dict(G.degree()) + node_orders = dict(G.degree()) from collections import defaultdict - counter = KeyedArray() + bond_order_counter = KeyedArray() for edge in G.edges(): - key = (min(degrees[edge[0]], degrees[edge[1]]), max(degrees[edge[0]], degrees[edge[1]])) - counter[key] += 1 + key = (min(node_orders[edge[0]], node_orders[edge[1]]), max(node_orders[edge[0]], node_orders[edge[1]])) + bond_order_counter[key] += 1 N = G.number_of_nodes() L = G.number_of_edges() + order_count = KeyedArray() + for pID, order in node_orders.items(): + order_count[order] += 1 + ## Write it to a file in output_dir #with open(configfilename + '_ChungLu.pkl', 'wb') as f: # pickle.dump({"N":N, "N_edges":L, "counters":counter}, f) # Calculate the modularity - return WeightedType(counter, 1) + retval = self.init() + retval["bond_order_count"] = WeightedType(bond_order_counter, 1) + retval["order_count"] = WeightedType(order_count, 1) + return retval def init(self): - return WeightedType(KeyedArray(), 0) - + return KeyedWeightedKeyedArray() +OutputFile.output_props["ChungLu"] = ChungLuConfigurationModel() -OutputFile.output_props["N"] = SingleAttrib('ParticleCount', 'val', [], [], [], missing_val=None)#We use missing_val=None to cause an error if the tag is missing -OutputFile.output_props["p"] = SingleAttrib('Pressure', 'Avg', [], [], [], missing_val=None) -OutputFile.output_props["cv"] = SingleAttrib('ResidualHeatCapacity', 'Value', [], [], [], div_by_N=True, missing_val=None) -OutputFile.output_props["u"] = SingleAttrib('UConfigurational', 'Mean', [], [], [], div_by_N=True, missing_val=None) -OutputFile.output_props["T"] = SingleAttrib('Temperature', 'Mean', [], [], [], missing_val=None) -OutputFile.output_props["density"] = SingleAttrib('Density', 'val', [], [], [], missing_val=None) -OutputFile.output_props["MSD"] = SingleAttrib('MSD/Species', 'diffusionCoeff', [], [], ['-LMSD'], missing_val=None, skip_missing=True) -OutputFile.output_props["NeventsSO"] = SingleAttrib('EventCounters/Entry[@Name="SOCells"]', # Outputfile tag name - 'Count', # Outputfile tag attribute name - ["Rso"], # Required state variable - [], # Required output variables - [], # Required output plugins - div_by_N=True, # Divide the count by N - div_by_t=True, # Also divide by t - missing_val=0) # If counter is missing, return 0 -OutputFile.output_props["VACF"] = VACFOutputProperty() -OutputFile.output_props["RadialDistEnd"] = RadialDistEndOutputProperty() -OutputFile.output_props["RadialDistribution"] = RadialDistributionOutputProperty() -OutputFile.output_props["FCCOrder"] = OrderParameterProperty(6) -OutputFile.output_props["CollisionMatrix"] = CollisionMatrixOutputProperty() -OutputFile.output_props["ChungLu"] = ChungLuConfigurationModel() \ No newline at end of file diff --git a/src/pydynamo/test_weighted_types.py b/src/pydynamo/test_weighted_types.py index b1ea9802d..8527b84ed 100644 --- a/src/pydynamo/test_weighted_types.py +++ b/src/pydynamo/test_weighted_types.py @@ -130,4 +130,23 @@ def test_keyed_keyed_array(): c = a * b assert c["a"]["1"] == pytest.approx(2) assert "b" not in c - assert "c" not in c \ No newline at end of file + assert "c" not in c + +def test_keyed_weighted_keyed_array(): + constructor = lambda : WeightedType(KeyedArray()) + constructor.__name__ = "WeightedType" + aval = WeightedType(KeyedArray(type=float, values={"1":1, "2": 2, "3":3}), 1) + bval = WeightedType(KeyedArray(type=float, values={"1":2, "2": 4}), 0.1) + a = KeyedArray(type=constructor, values = { + "a": aval, + }) + + b = KeyedArray(type=constructor, values = { + "a": bval, + }) + + c = a + b + + assert c["a"].avg()["1"] == pytest.approx((1 * 1 + 2 * 0.1) / (1 + 0.1)) + assert c["a"].avg()["2"] == pytest.approx((2 * 1 + 4 * 0.1) / (1 + 0.1)) + assert c["a"].avg()["3"] == pytest.approx((3 * 1 + 0 * 0.1) / (1 + 0.1)) \ No newline at end of file diff --git a/src/pydynamo/weighted_types.py b/src/pydynamo/weighted_types.py index 392da5806..3dd0301a1 100644 --- a/src/pydynamo/weighted_types.py +++ b/src/pydynamo/weighted_types.py @@ -149,7 +149,7 @@ class WeightedType(): '''This class implements weighted arithmetic means along with an estimate of the standard error for the mean. ''' - def __init__(self, value = 0, weight = 0): + def __init__(self, value = 0.0, weight = 0.0): """Initialise the weighted value.""" ww = weight * weight vv = element_wise_multiply(value, value) @@ -163,10 +163,11 @@ def __init__(self, value = 0, weight = 0): def __add__(self, v): if not isinstance(v, WeightedType): - if not isinstance(v._w_v_sum, type(self._w_v_sum)): - raise RuntimeError("Cannot add non-WeightedType to WeightedType") raise RuntimeError("Cannot add non-WeightedType to WeightedType") + if not isinstance(v._w_v_sum, type(self._w_v_sum)): + raise RuntimeError(f"WeightedTypes have incompatible types {type(self._w_v_sum)} and {type(v._w_v_sum)}") + import copy retval = copy.copy(v) if self._w_sum == 0: @@ -223,7 +224,7 @@ def ufloat(self): # If the average is an array, we need to convert it to a ufloat array return uncertainties.unumpy.uarray(avg, std_dev) elif isinstance(avg, KeyedArray): - return KeyedArray({k: uncertainties.ufloat(avg[k], std_dev[k]) for k in avg.keys()}) + return KeyedArray(uncertainties.ufloat, values={k: uncertainties.ufloat(avg[k], std_dev[k]) for k in avg.keys()}) else: return uncertainties.ufloat(avg, std_dev) @@ -233,4 +234,15 @@ def __str__(self): def __repr__(self): return repr(self.ufloat()) + +class WeightedKeyedArray(WeightedType): + """This type is picklable and can be used in multiprocessing. + """ + def __init__(self, type = float, values = None): + super().__init__(KeyedArray()) + +class KeyedWeightedKeyedArray(KeyedArray): + def __init__(self): + super().__init__(type=WeightedKeyedArray) + from collections import defaultdict From f6611a0d62f7cc71262ca4fc35dcabd6cbd4c9e9 Mon Sep 17 00:00:00 2001 From: Marcus Bannerman Date: Thu, 1 May 2025 13:46:37 +0100 Subject: [PATCH 08/19] Changes for plots for leo --- scripts/plotter.ipynb | 1 + 1 file changed, 1 insertion(+) diff --git a/scripts/plotter.ipynb b/scripts/plotter.ipynb index 338ff45d9..f590192ff 100644 --- a/scripts/plotter.ipynb +++ b/scripts/plotter.ipynb @@ -69,6 +69,7 @@ " )\n", "plt.xlabel(\"$k_i\\\\,k_j$\")\n", "plt.ylabel(\"$\\\\frac{N(k_i,\\\\,k_j)}{N(k_i)\\\\,N(k_j)}$\")\n", + "plt.yscale(\"log\")\n", "plt.legend()\n", "plt.grid()\n", "plt.show()" From 0cde04796f172df6ce3f7fda3119ae09a6bfc76a Mon Sep 17 00:00:00 2001 From: Marcus Bannerman Date: Thu, 5 Jun 2025 15:20:21 +0100 Subject: [PATCH 09/19] Typo fix --- src/dynamo/dynamo/outputplugins/collMatrix.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/dynamo/dynamo/outputplugins/collMatrix.hpp b/src/dynamo/dynamo/outputplugins/collMatrix.hpp index c967e4529..c4e877806 100644 --- a/src/dynamo/dynamo/outputplugins/collMatrix.hpp +++ b/src/dynamo/dynamo/outputplugins/collMatrix.hpp @@ -51,7 +51,7 @@ class OPCollMatrix : public OutputPlugin { unsigned long totalCount; - // EventKet is a pair of EventSourceKey and EEventType + // EventKey is a pair of EventSourceKey and EEventType // It describes the type of event and its source //! A key for two events From 8a5afd3066b67b8d83c2f2fa1bbafdf75972c7a2 Mon Sep 17 00:00:00 2001 From: Marcus Bannerman Date: Fri, 25 Jul 2025 13:21:20 +0100 Subject: [PATCH 10/19] Added multiple outputs from outputplugins to be processed --- .gitignore | 4 +- src/pydynamo/__init__.py | 21 ++++++--- src/pydynamo/output_properties.py | 75 +++++++++++++------------------ 3 files changed, 49 insertions(+), 51 deletions(-) diff --git a/.gitignore b/.gitignore index 9f98cc17f..120db7184 100644 --- a/.gitignore +++ b/.gitignore @@ -12,4 +12,6 @@ __pycache__ Testing/Temporary/CTestCostData.txt .eggs wheelhouse -vcpkg_installed \ No newline at end of file +vcpkg_installed + +result diff --git a/src/pydynamo/__init__.py b/src/pydynamo/__init__.py index 13a261957..ba7429e68 100755 --- a/src/pydynamo/__init__.py +++ b/src/pydynamo/__init__.py @@ -185,12 +185,21 @@ def fetch_data_worker(args): dataout["tTotal"] += outputfile.t() for prop in manager.outputs: - outputplugin = OutputFile.output_props[prop] - result = outputplugin.result(state, outputfile, configfilename, counter, manager, output_dir) - if result != None: - if prop not in dataout: - dataout[prop] = outputplugin.init() - dataout[prop] += result + try: + outputplugin = OutputFile.output_props[prop] + # Output plugins will return dictionaries of properties to let them return multiple properties + + result = outputplugin.result(state, outputfile, configfilename, counter, manager, output_dir) + initValues = outputplugin.init() + for propname in result: + # Ensure its initialised if needed + if propname not in dataout: + dataout[propname] = initValues[propname] + dataout[propname] += result[propname] + except Exception as e: + print("Error while processing output property", prop, "in", output_dir, ":", e) + raise + #except Exception as e: # print("Processing", output_dir, " gave exception", e) # #raise diff --git a/src/pydynamo/output_properties.py b/src/pydynamo/output_properties.py index dec71db1d..1a5a91f9f 100644 --- a/src/pydynamo/output_properties.py +++ b/src/pydynamo/output_properties.py @@ -49,14 +49,15 @@ def __init__(self, dependent_statevars : list, dependent_outputs : list, depende self._dep_outputplugins = dependent_outputplugins def init(self): - return None + return {} def result(self, state, outputfile, configfilename, counter, manager, output_dir): - return None + return {} class SingleAttrib(OutputProperty): - def __init__(self, tag, attrib, dependent_statevars, dependent_outputs, dependent_outputplugins, time_weighted=True, div_by_N=False, div_by_t=False, missing_val = 0, skip_missing=False): + def __init__(self, propkey, tag, attrib, dependent_statevars, dependent_outputs, dependent_outputplugins, time_weighted=True, div_by_N=False, div_by_t=False, missing_val = 0, skip_missing=False): OutputProperty.__init__(self, dependent_statevars, dependent_outputs, dependent_outputplugins) + self._propkey = propkey self._tag = tag self._attrib = attrib self._time_weighted = time_weighted @@ -66,7 +67,7 @@ def __init__(self, tag, attrib, dependent_statevars, dependent_outputs, dependen self._skip_missing=skip_missing def init(self): - return WeightedType() + return {self._propkey: WeightedType()} def value(self, outputfile): tag = outputfile.tree.find('.//'+self._tag) @@ -105,16 +106,16 @@ def weight(self, outputfile): return float(outputfile.tree.find('.//Duration').attrib['Events']) def result(self, state, outputfile, configfilename, counter, manager, output_dir): - return WeightedType(self.value(outputfile), self.weight(outputfile)) - -OutputFile.output_props["N"] = SingleAttrib('ParticleCount', 'val', [], [], [], missing_val=None)#We use missing_val=None to cause an error if the tag is missing -OutputFile.output_props["p"] = SingleAttrib('Pressure', 'Avg', [], [], [], missing_val=None) -OutputFile.output_props["cv"] = SingleAttrib('ResidualHeatCapacity', 'Value', [], [], [], div_by_N=True, missing_val=None) -OutputFile.output_props["u"] = SingleAttrib('UConfigurational', 'Mean', [], [], [], div_by_N=True, missing_val=None) -OutputFile.output_props["T"] = SingleAttrib('Temperature', 'Mean', [], [], [], missing_val=None) -OutputFile.output_props["density"] = SingleAttrib('Density', 'val', [], [], [], missing_val=None) -OutputFile.output_props["MSD"] = SingleAttrib('MSD/Species', 'diffusionCoeff', [], [], ['-LMSD'], missing_val=None, skip_missing=True) -OutputFile.output_props["NeventsSO"] = SingleAttrib('EventCounters/Entry[@Name="SOCells"]', # Outputfile tag name + return {self._propkey: WeightedType(self.value(outputfile), self.weight(outputfile))} + +OutputFile.output_props["N"] = SingleAttrib("N", 'ParticleCount', 'val', [], [], [], missing_val=None)#We use missing_val=None to cause an error if the tag is missing +OutputFile.output_props["p"] = SingleAttrib("p", 'Pressure', 'Avg', [], [], [], missing_val=None) +OutputFile.output_props["cv"] = SingleAttrib("cv", 'ResidualHeatCapacity', 'Value', [], [], [], div_by_N=True, missing_val=None) +OutputFile.output_props["u"] = SingleAttrib("u", 'UConfigurational', 'Mean', [], [], [], div_by_N=True, missing_val=None) +OutputFile.output_props["T"] = SingleAttrib("T",'Temperature', 'Mean', [], [], [], missing_val=None) +OutputFile.output_props["density"] = SingleAttrib("density", 'Density', 'val', [], [], [], missing_val=None) +OutputFile.output_props["MSD"] = SingleAttrib("MSD",'MSD/Species', 'diffusionCoeff', [], [], ['-LMSD'], missing_val=None, skip_missing=True) +OutputFile.output_props["NeventsSO"] = SingleAttrib("NeventsSO", 'EventCounters/Entry[@Name="SOCells"]', # Outputfile tag name 'Count', # Outputfile tag attribute name ["Rso"], # Required state variable [], # Required output variables @@ -134,9 +135,6 @@ def parseToArray(text): class CollisionMatrixOutputProperty(OutputProperty): def __init__(self): OutputProperty.__init__(self, dependent_statevars=[], dependent_outputs=[], dependent_outputplugins=['-LCollisionMatrix']) - - def result(self, state, outputfile, configfilename, counter, manager, output_dir): - return None OutputFile.output_props["CollisionMatrix"] = CollisionMatrixOutputProperty() class VACFOutputProperty(OutputProperty): @@ -152,15 +150,16 @@ def result(self, state, outputfile, configfilename, counter, manager, output_dir for tag in outputfile.tree.findall('.//VACF/Topology/Structure'): pickle.dump(parseToArray(tag.text), open(filename_root + '/topology_'+tag.attrib['Name']+'.pkl', 'wb')) - return None + return {} OutputFile.output_props["VACF"] = VACFOutputProperty() + class RadialDistributionOutputProperty(OutputProperty): def __init__(self): OutputProperty.__init__(self, dependent_statevars=[], dependent_outputs=[], dependent_outputplugins=['-LRadialDistribution']) def init(self): - return WeightedType() + return {"RadialDistribution":WeightedType()} def result(self, state, outputfile, configfilename, counter, manager, output_dir): #Presume that each tag is in order, and has a common bin width @@ -200,7 +199,7 @@ def result(self, state, outputfile, configfilename, counter, manager, output_dir for n in range(2, moments.shape[0]): central_moments[n-1] = sum(scipy.special.comb(n, i) * moments[i] * (N0-central_moments[0,:])**(n-i) for i in range(n+1)) - return WeightedType(central_moments, samples) + return {"RadialDistribution":WeightedType(central_moments, samples)} OutputFile.output_props["RadialDistribution"] = RadialDistributionOutputProperty() class RadialDistEndOutputProperty(OutputProperty): @@ -215,7 +214,7 @@ def result(self, state, outputfile, configfilename, counter, manager, output_dir filename_root = manager.workdir+'_RadialDist/' + manager.statename(state, var_separator='/') if os.path.isdir(filename_root): - return None + return {} os.makedirs(filename_root, exist_ok=True) @@ -228,7 +227,7 @@ def result(self, state, outputfile, configfilename, counter, manager, output_dir B = tag.attrib['Name2'] output_pkl = filename_root + '/species_'+A+'_'+B+'.pkl' pickle.dump(parseToArray(tag.text), open(output_pkl, 'wb')) - return None + return {} OutputFile.output_props["RadialDistEnd"] = RadialDistEndOutputProperty() class OrderParameterProperty(OutputProperty): @@ -238,9 +237,8 @@ class OrderParameterProperty(OutputProperty): This property is expensive to run at data collection time, as it processes every configuration file to determine the order parameter. The advantage is that it can be run on any simulation, no need for extra output plugins. ''' - def __init__(self, L): + def __init__(self): OutputProperty.__init__(self, dependent_statevars=[], dependent_outputs=[], dependent_outputplugins=[]) - self.L = L def result(self, state, outputfile, configfilename, counter, manager, output_dir): # We use freud to calculate the Steinhardt order parameter @@ -251,27 +249,16 @@ def result(self, state, outputfile, configfilename, counter, manager, output_dir box, points = configfile.to_freud() #Steinhardt for FCC - ql = freud.order.Steinhardt(self.L) - ql.compute((box, points), neighbors={"num_neighbors": self.L}) + L=6 + ql = freud.order.Steinhardt(L) + ql.compute((box, points), neighbors={"num_neighbors": L}) ql_value = ql.particle_order - return WeightedType(numpy.mean(ql_value), 1) + return {"FCCOrder":WeightedType(numpy.mean(ql_value), 1)} def init(self): - return WeightedType() - - def result(self, state, outputfile, configfilename, counter, manager, output_dir): - import freud - configfile = ConfigFile(configfilename) - box, points = configfile.to_freud() - - #Steinhardt for FCC - ql = freud.order.Steinhardt(self.L) - ql.compute((box, points), neighbors={"num_neighbors": self.L}) - ql_value = ql.particle_order - - return WeightedType(numpy.mean(ql_value), 1) -OutputFile.output_props["FCCOrder"] = OrderParameterProperty(6) + return {"FCCOrder":WeightedType()} +OutputFile.output_props["OrderParameter"] = OrderParameterProperty() class ChungLuConfigurationModel(OutputProperty): ''' @@ -319,11 +306,11 @@ def result(self, state, outputfile, configfilename, bond_order_counter, manager, # Calculate the modularity retval = self.init() - retval["bond_order_count"] = WeightedType(bond_order_counter, 1) - retval["order_count"] = WeightedType(order_count, 1) + retval["ChungLu"]["bond_order_count"] = WeightedType(bond_order_counter, 1) + retval["ChungLu"]["order_count"] = WeightedType(order_count, 1) return retval def init(self): - return KeyedWeightedKeyedArray() + return {"ChungLu":KeyedWeightedKeyedArray()} OutputFile.output_props["ChungLu"] = ChungLuConfigurationModel() From b137519d50e04eb5c7ffd516b26bea76d54e32de Mon Sep 17 00:00:00 2001 From: Marcus Bannerman Date: Fri, 25 Jul 2025 14:40:22 +0100 Subject: [PATCH 11/19] Pulled in final data types for collision matrix --- src/pydynamo/output_properties.py | 36 +++++++++++++- src/pydynamo/weighted_types.py | 83 +++++++++++++++++++++++++++++++ 2 files changed, 117 insertions(+), 2 deletions(-) diff --git a/src/pydynamo/output_properties.py b/src/pydynamo/output_properties.py index 1a5a91f9f..c15367d5a 100644 --- a/src/pydynamo/output_properties.py +++ b/src/pydynamo/output_properties.py @@ -7,8 +7,8 @@ from pydynamo.config_files import ConfigFile from pydynamo.file_types import XMLFile, validate_xmlfile -from pydynamo.weighted_types import (KeyedArray, KeyedWeightedKeyedArray, - WeightedType) +from pydynamo.weighted_types import (KeyedArray, KeyedWeightedKeyedArray, KeyedKeyedArray, + WeightedType, Histogram) # A XMLFile/ElementTree but specialised for DynamO output files @@ -314,3 +314,35 @@ def init(self): return {"ChungLu":KeyedWeightedKeyedArray()} OutputFile.output_props["ChungLu"] = ChungLuConfigurationModel() +class CollisionMatrixOutputProperty(OutputProperty): + def __init__(self): + OutputProperty.__init__(self, dependent_statevars=[], dependent_outputs=[], dependent_outputplugins=['-LCollisionMatrix']) + + def result(self, state, outputfile: OutputFile, configfilename, counter, manager, output_dir): + N = outputfile.N() + t = outputfile.t() + + retval = self.init() + retval["CaptureStateHistogram"] = Histogram.load(outputfile.tree.find(".//CaptureStateHistogram/HistogramWeighted")) + + Rates = KeyedKeyedArray() + for tag in outputfile.tree.findall(".//CollCounters/PairCaptureCounters/Count"): + iName = tag.attrib["Name"] + eType = tag.attrib["Event"] + minCap = int(tag.attrib["captures1"]) + maxCap = int(tag.attrib["captures2"]) + count = int(tag.find("./RijDotVij/Histogram").attrib["SampleCount"]) / 2 + Rates[iName+"_"+eType][(minCap, maxCap)] = count / N / t + + for key, item in Rates.items(): + retval["CollisionMatrix"][key] = WeightedType(item, outputfile.t()) + + return retval + + def init(self): + return { + "CollisionMatrix":KeyedWeightedKeyedArray(), + "CaptureStateHistogram": Histogram() + } + +OutputFile.output_props["CollisionMatrix"] = CollisionMatrixOutputProperty() diff --git a/src/pydynamo/weighted_types.py b/src/pydynamo/weighted_types.py index 3dd0301a1..3af924e44 100644 --- a/src/pydynamo/weighted_types.py +++ b/src/pydynamo/weighted_types.py @@ -245,4 +245,87 @@ class KeyedWeightedKeyedArray(KeyedArray): def __init__(self): super().__init__(type=WeightedKeyedArray) +class KeyedKeyedArray(KeyedArray): + def __init__(self): + super().__init__(type=KeyedArray) + +class Histogram: + """ + A histogram class to hold the histogram data. + """ + def __init__(self, binwidth = None, total_weight = None): + self.binwidth = binwidth + self.total_weight = total_weight + self.data = KeyedArray() + + def insert(self, key, value): + self.data[key] = value + + def __getitem__(self, key): + return self.data[key] + + def __repr__(self): + return f"Histogram(binwidth={self.binwidth}, total_weight={self.total_weight}, data={self.data})" + + def __add__(self, other): + if not isinstance(other, Histogram): + raise TypeError("Can only add another Histogram") + + if self.binwidth is None: + return other + if other.binwidth is None: + return self + if self.binwidth != other.binwidth: + raise ValueError("Cannot add Histograms with different bin widths") + + new_histogram = Histogram(self.binwidth, self.total_weight + other.total_weight) + for key in set(self.keys()).union(other.keys()): + new_histogram.insert(key, self.get(key, 0) + other.get(key, 0)) + + return new_histogram + + def get(self, key, default=None): + """ + Get the value for a key, or default if not found. + """ + if key in self.data: + return self.data[key] + if default is None: + raise KeyError(f"Key {key} not found in histogram") + return default + + def keys(self): + """ + Return the keys of the histogram. + """ + return self.data.keys() + + @staticmethod + def load(tag): + """ + Load a histogram from an XML tag. + """ + if tag is None: + raise ValueError("Histogram is None") + + if "TotalWeight" in tag.attrib: + weight = float(tag.attrib["TotalWeight"]) + elif "SampleCount" in tag.attrib: + weight = float(tag.attrib["SampleCount"]) + else: + raise ValueError("Histogram tag must have 'TotalWeight' or 'SampleCount' attributes.") + + binwidth = float(tag.attrib["BinWidth"]) + + retval = Histogram(binwidth, weight) + + if tag.text is not None: + for line in tag.text.splitlines(): + if line.strip() == "": + continue + key, value= list(map(float, line.strip().split())) + retval.insert(key, value) + + return retval + from collections import defaultdict From 787d7768fefe23fa92544621f42c00205f5eac93 Mon Sep 17 00:00:00 2001 From: Marcus Bannerman Date: Fri, 25 Jul 2025 18:57:36 +0100 Subject: [PATCH 12/19] Fix for histograms to load and keep weight --- src/pydynamo/weighted_types.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/pydynamo/weighted_types.py b/src/pydynamo/weighted_types.py index 3af924e44..894025cc3 100644 --- a/src/pydynamo/weighted_types.py +++ b/src/pydynamo/weighted_types.py @@ -324,7 +324,7 @@ def load(tag): if line.strip() == "": continue key, value= list(map(float, line.strip().split())) - retval.insert(key, value) + retval.insert(key, value * weight) return retval From 6c59a2eb5d553385ae5a4a1e6946511828a3f781 Mon Sep 17 00:00:00 2001 From: Marcus Bannerman Date: Sat, 9 Aug 2025 14:38:24 +0100 Subject: [PATCH 13/19] Made the visualiser optional in nix --- derivation.nix | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/derivation.nix b/derivation.nix index 1ff59a038..f92d0df2f 100644 --- a/derivation.nix +++ b/derivation.nix @@ -2,7 +2,9 @@ # To install `nix-env -u -f default.nix` # To develop `nix-shell` (will build the shell with dependencies) # To test build `nix-build` -{ pkgs, python3 }: +{ pkgs, python3, + visualiser ? true, + }: python3.pkgs.buildPythonPackage rec { name = "pydynamo"; src = ./.; @@ -29,15 +31,17 @@ python3.pkgs.buildPythonPackage rec { gcc pkg-config clang-tools + ] ++ propagatedBuildInputs + ++ (lib.optionals visualiser [ wrapGAppsHook3 - ] ++ propagatedBuildInputs; + ]); buildInputs = with pkgs; [ # Basic build dependencies bzip2.dev boost.dev eigen - + ] ++ (lib.optionals visualiser [ # Visualiser libGL gtkmm3.dev @@ -47,5 +51,5 @@ python3.pkgs.buildPythonPackage rec { cairomm.dev libpng mesa - ]; + ]); } From 5101238569916b4dfbe54e4baf106612cf5ce336 Mon Sep 17 00:00:00 2001 From: Marcus Bannerman Date: Mon, 11 Aug 2025 00:26:11 +0100 Subject: [PATCH 14/19] Added the V2 output to chunglu --- src/pydynamo/output_properties.py | 19 +++++++++++++++---- src/pydynamo/weighted_types.py | 24 +++++++++++++++++++----- 2 files changed, 34 insertions(+), 9 deletions(-) diff --git a/src/pydynamo/output_properties.py b/src/pydynamo/output_properties.py index c15367d5a..53ac838ce 100644 --- a/src/pydynamo/output_properties.py +++ b/src/pydynamo/output_properties.py @@ -8,7 +8,7 @@ from pydynamo.config_files import ConfigFile from pydynamo.file_types import XMLFile, validate_xmlfile from pydynamo.weighted_types import (KeyedArray, KeyedWeightedKeyedArray, KeyedKeyedArray, - WeightedType, Histogram) + WeightedType, Histogram, KeyedHistogram) # A XMLFile/ElementTree but specialised for DynamO output files @@ -323,8 +323,11 @@ def result(self, state, outputfile: OutputFile, configfilename, counter, manager t = outputfile.t() retval = self.init() + + # First, parse out the histogram of capture states retval["CaptureStateHistogram"] = Histogram.load(outputfile.tree.find(".//CaptureStateHistogram/HistogramWeighted")) + # Next, get the per-particle rates of events between various capture levels and event types. Rates = KeyedKeyedArray() for tag in outputfile.tree.findall(".//CollCounters/PairCaptureCounters/Count"): iName = tag.attrib["Name"] @@ -333,16 +336,24 @@ def result(self, state, outputfile: OutputFile, configfilename, counter, manager maxCap = int(tag.attrib["captures2"]) count = int(tag.find("./RijDotVij/Histogram").attrib["SampleCount"]) / 2 Rates[iName+"_"+eType][(minCap, maxCap)] = count / N / t - for key, item in Rates.items(): retval["CollisionMatrix"][key] = WeightedType(item, outputfile.t()) + # Now get the properties for a particular capture count + for tag in outputfile.tree.findall(".//CollCounters/CaptureCounters/Count"): + iName = tag.attrib["Name"] + eType = tag.attrib["Event"] + captures = int(tag.attrib["captures"]) + # Square velocity + retval["V2"][(iName, eType, captures)] += Histogram.load(tag.find("./V2/Histogram")) + return retval def init(self): return { - "CollisionMatrix":KeyedWeightedKeyedArray(), - "CaptureStateHistogram": Histogram() + "CollisionMatrix": KeyedWeightedKeyedArray(), + "CaptureStateHistogram": Histogram(), + "V2": KeyedHistogram(), } OutputFile.output_props["CollisionMatrix"] = CollisionMatrixOutputProperty() diff --git a/src/pydynamo/weighted_types.py b/src/pydynamo/weighted_types.py index 894025cc3..546739915 100644 --- a/src/pydynamo/weighted_types.py +++ b/src/pydynamo/weighted_types.py @@ -253,9 +253,8 @@ class Histogram: """ A histogram class to hold the histogram data. """ - def __init__(self, binwidth = None, total_weight = None): + def __init__(self, binwidth = None): self.binwidth = binwidth - self.total_weight = total_weight self.data = KeyedArray() def insert(self, key, value): @@ -265,7 +264,7 @@ def __getitem__(self, key): return self.data[key] def __repr__(self): - return f"Histogram(binwidth={self.binwidth}, total_weight={self.total_weight}, data={self.data})" + return f"Histogram(binwidth={self.binwidth}, data={self.data})" def __add__(self, other): if not isinstance(other, Histogram): @@ -278,11 +277,17 @@ def __add__(self, other): if self.binwidth != other.binwidth: raise ValueError("Cannot add Histograms with different bin widths") - new_histogram = Histogram(self.binwidth, self.total_weight + other.total_weight) + new_histogram = Histogram(self.binwidth) for key in set(self.keys()).union(other.keys()): new_histogram.insert(key, self.get(key, 0) + other.get(key, 0)) return new_histogram + + def total_weight(self): + return sum([value for _, value in self.data.items()]) * self.binwidth + + def avg(): + return sum([key * value for key, value in self.data.items()]) / sum([value for _, value in self.data.items()]) def get(self, key, default=None): """ @@ -317,15 +322,24 @@ def load(tag): binwidth = float(tag.attrib["BinWidth"]) - retval = Histogram(binwidth, weight) + retval = Histogram(binwidth) + value_sum = 0 if tag.text is not None: for line in tag.text.splitlines(): if line.strip() == "": continue key, value= list(map(float, line.strip().split())) retval.insert(key, value * weight) + value_sum += value + if (value_sum > 0) and abs(value_sum * retval.binwidth - 1.0) > 1e-3: + print("Warning! Unnormalised histogram", value_sum * retval.binwidth) + return retval +class KeyedHistogram(KeyedArray): + def __init__(self): + super().__init__(type=Histogram) + from collections import defaultdict From 1f9c2b01ec502f8011ffb64dbf7006146bb1700f Mon Sep 17 00:00:00 2001 From: Marcus Bannerman Date: Mon, 11 Aug 2025 21:42:10 +0100 Subject: [PATCH 15/19] Made histograms used WeightedKeyedArrays to preserve uncertainty and added a correct average calculator --- src/pydynamo/weighted_types.py | 40 ++++++++++++++++------------------ 1 file changed, 19 insertions(+), 21 deletions(-) diff --git a/src/pydynamo/weighted_types.py b/src/pydynamo/weighted_types.py index 546739915..a7338c049 100644 --- a/src/pydynamo/weighted_types.py +++ b/src/pydynamo/weighted_types.py @@ -238,8 +238,8 @@ def __repr__(self): class WeightedKeyedArray(WeightedType): """This type is picklable and can be used in multiprocessing. """ - def __init__(self, type = float, values = None): - super().__init__(KeyedArray()) + def __init__(self, type = float, values = None, weight = 0.0): + super().__init__(value= KeyedArray() if (values == None) else values, weight=weight) class KeyedWeightedKeyedArray(KeyedArray): def __init__(self): @@ -255,13 +255,7 @@ class Histogram: """ def __init__(self, binwidth = None): self.binwidth = binwidth - self.data = KeyedArray() - - def insert(self, key, value): - self.data[key] = value - - def __getitem__(self, key): - return self.data[key] + self.data = WeightedKeyedArray() def __repr__(self): return f"Histogram(binwidth={self.binwidth}, data={self.data})" @@ -278,17 +272,19 @@ def __add__(self, other): raise ValueError("Cannot add Histograms with different bin widths") new_histogram = Histogram(self.binwidth) - for key in set(self.keys()).union(other.keys()): - new_histogram.insert(key, self.get(key, 0) + other.get(key, 0)) - + new_histogram.data = self.data + other.data return new_histogram - - def total_weight(self): - return sum([value for _, value in self.data.items()]) * self.binwidth - - def avg(): - return sum([key * value for key, value in self.data.items()]) / sum([value for _, value in self.data.items()]) + def avg(self): + udata = self.data.ufloat() + sum = 0.0 + vsum = 0.0 + for k, v in udata.items(): + sum += k * v + vsum += v + assert abs(vsum * self.binwidth - 1.0) < 1e-3, "Histogram is not normalised!" + return sum * self.binwidth + def get(self, key, default=None): """ Get the value for a key, or default if not found. @@ -322,17 +318,19 @@ def load(tag): binwidth = float(tag.attrib["BinWidth"]) - retval = Histogram(binwidth) - + values = KeyedArray() value_sum = 0 if tag.text is not None: for line in tag.text.splitlines(): if line.strip() == "": continue key, value= list(map(float, line.strip().split())) - retval.insert(key, value * weight) + values[key] = value value_sum += value + retval = Histogram(binwidth) + retval.data = WeightedKeyedArray(type=float, values=values, weight=weight) + if (value_sum > 0) and abs(value_sum * retval.binwidth - 1.0) > 1e-3: print("Warning! Unnormalised histogram", value_sum * retval.binwidth) From aa276e8db6cccb27340c6e48f58e67cd0c73c071 Mon Sep 17 00:00:00 2001 From: Marcus Bannerman Date: Wed, 13 Aug 2025 12:59:52 +0100 Subject: [PATCH 16/19] Added an optimised histogram average --- src/pydynamo/weighted_types.py | 26 ++++++++++++++++++-------- 1 file changed, 18 insertions(+), 8 deletions(-) diff --git a/src/pydynamo/weighted_types.py b/src/pydynamo/weighted_types.py index a7338c049..5af2ad81b 100644 --- a/src/pydynamo/weighted_types.py +++ b/src/pydynamo/weighted_types.py @@ -186,6 +186,7 @@ def __add__(self, v): return retval def stats(self): + # Note, there's an implementation for keyed array averages in the Histogram class. if self._w_sum == 0: if isinstance(self._w_v_sum, numpy.ndarray): v = numpy.empty_like(self._w_v_sum) @@ -276,14 +277,23 @@ def __add__(self, other): return new_histogram def avg(self): - udata = self.data.ufloat() - sum = 0.0 - vsum = 0.0 - for k, v in udata.items(): - sum += k * v - vsum += v - assert abs(vsum * self.binwidth - 1.0) < 1e-3, "Histogram is not normalised!" - return sum * self.binwidth + # This is a slow method using uncertainties, but it is correct. + #udata = self.data.ufloat() + #sum = 0.0 + #vsum = 0.0 + #for k, v in udata.items(): + # sum += k * v + # vsum += v + #assert abs(vsum * self.binwidth - 1.0) < 1e-3, "Histogram is not normalised!" + #return sum * self.binwidth + + # About 30x faster than the above method + vavg = self.data._w_v_sum / self.data._w_sum + avg = sum([k * v for k,v in vavg.items()]) * self.binwidth + vavg_var = (1 / self.data._count) * (self.data._w_sum * element_wise_multiply(vavg, vavg) - 2 * element_wise_multiply(vavg, self.data._w_v_sum) + self.data._w_vv_sum) / self.data._w_sum + avg_std_dev = math.sqrt(sum([k * k * v for k, v in vavg_var.items()])) * self.binwidth + return uncertainties.ufloat(avg, avg_std_dev) + def get(self, key, default=None): """ From cd3c7c5a63fd3b5152723af477277533d0fc7b8f Mon Sep 17 00:00:00 2001 From: Marcus Bannerman Date: Fri, 4 Sep 2026 12:54:23 +0100 Subject: [PATCH 17/19] Added system MFT histograms --- .gitignore | 3 + CMakeLists.txt | 12 +- derivation.nix | 6 +- flake.lock | 61 +++++++ flake.nix | 26 +++ scripts/HS_stats.py | 155 ++++++++++++++++++ .../dynamo/outputplugins/collMatrix.cpp | 44 ++++- .../dynamo/outputplugins/collMatrix.hpp | 15 +- .../outputplugins/eventtypetracking.hpp | 4 +- 9 files changed, 312 insertions(+), 14 deletions(-) create mode 100644 flake.lock create mode 100644 flake.nix create mode 100755 scripts/HS_stats.py diff --git a/.gitignore b/.gitignore index 120db7184..a58805ca4 100644 --- a/.gitignore +++ b/.gitignore @@ -15,3 +15,6 @@ wheelhouse vcpkg_installed result +build_* +result +result-* \ No newline at end of file diff --git a/CMakeLists.txt b/CMakeLists.txt index eacc8ab17..44543635a 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -101,14 +101,18 @@ endif() ######### Dependencies ###################################################################### -set(Boost_USE_STATIC_LIBS ON) set(Boost_USE_MULTITHREADED ON) set(Boost_USE_STATIC_RUNTIME OFF) -find_package(Boost 1.66.0 COMPONENTS system filesystem program_options unit_test_framework) + +# Try static Boost first, then fallback to shared Boost +find_package(Boost 1.66.0 COMPONENTS filesystem program_options unit_test_framework) if(NOT Boost_FOUND) - message(WARNING "Static Boost not found, trying shared Boost") set(Boost_USE_STATIC_LIBS OFF) - find_package(Boost 1.66.0 REQUIRED COMPONENTS system filesystem program_options unit_test_framework) + find_package(Boost 1.66.0 REQUIRED COMPONENTS filesystem program_options unit_test_framework) +endif() + +if(NOT TARGET Boost::system) + add_library(Boost::system INTERFACE IMPORTED) endif() add_subdirectory(src/magnet) diff --git a/derivation.nix b/derivation.nix index f92d0df2f..90d802990 100644 --- a/derivation.nix +++ b/derivation.nix @@ -3,7 +3,7 @@ # To develop `nix-shell` (will build the shell with dependencies) # To test build `nix-build` { pkgs, python3, - visualiser ? true, + visualiser ? false, }: python3.pkgs.buildPythonPackage rec { name = "pydynamo"; @@ -38,8 +38,8 @@ python3.pkgs.buildPythonPackage rec { buildInputs = with pkgs; [ # Basic build dependencies - bzip2.dev - boost.dev + bzip2 + boost eigen ] ++ (lib.optionals visualiser [ # Visualiser diff --git a/flake.lock b/flake.lock new file mode 100644 index 000000000..032fb8b8e --- /dev/null +++ b/flake.lock @@ -0,0 +1,61 @@ +{ + "nodes": { + "flake-utils": { + "inputs": { + "systems": "systems" + }, + "locked": { + "lastModified": 1731533236, + "narHash": "sha256-l0KFg5HjrsfsO/JpG+r7fRrqm12kzFHyUHqHCVpMMbI=", + "owner": "numtide", + "repo": "flake-utils", + "rev": "11707dc2f618dd54ca8739b309ec4fc024de578b", + "type": "github" + }, + "original": { + "owner": "numtide", + "repo": "flake-utils", + "type": "github" + } + }, + "nixpkgs": { + "locked": { + "lastModified": 1751274312, + "narHash": "sha256-/bVBlRpECLVzjV19t5KMdMFWSwKLtb5RyXdjz3LJT+g=", + "owner": "NixOS", + "repo": "nixpkgs", + "rev": "50ab793786d9de88ee30ec4e4c24fb4236fc2674", + "type": "github" + }, + "original": { + "owner": "NixOS", + "ref": "nixos-24.11", + "repo": "nixpkgs", + "type": "github" + } + }, + "root": { + "inputs": { + "flake-utils": "flake-utils", + "nixpkgs": "nixpkgs" + } + }, + "systems": { + "locked": { + "lastModified": 1681028828, + "narHash": "sha256-Vy1rq5AaRuLzOxct8nz4T6wlgyUR7zLU309k9mBC768=", + "owner": "nix-systems", + "repo": "default", + "rev": "da67096a3b9bf56a91d16901293e51ba5b49a27e", + "type": "github" + }, + "original": { + "owner": "nix-systems", + "repo": "default", + "type": "github" + } + } + }, + "root": "root", + "version": 7 +} diff --git a/flake.nix b/flake.nix new file mode 100644 index 000000000..8207f305f --- /dev/null +++ b/flake.nix @@ -0,0 +1,26 @@ +{ + description = "PyDynamO / DynamO build environment"; + + inputs = { + nixpkgs.url = "github:NixOS/nixpkgs/nixos-24.11"; + flake-utils.url = "github:numtide/flake-utils"; + }; + + outputs = { self, nixpkgs, flake-utils }: + flake-utils.lib.eachDefaultSystem (system: + let + pkgs = import nixpkgs { + inherit system; + }; + pydynamo = pkgs.callPackage ./derivation.nix {}; + in + { + packages.default = pydynamo; + packages.pydynamo = pydynamo; + + devShells.default = pkgs.mkShell { + inputsFrom = [ pydynamo ]; + }; + } + ); +} diff --git a/scripts/HS_stats.py b/scripts/HS_stats.py new file mode 100755 index 000000000..6591ebb08 --- /dev/null +++ b/scripts/HS_stats.py @@ -0,0 +1,155 @@ +#!/usr/bin/env python3 +import pydynamo +from pydynamo import ET + + +def setup_worker( config, # The name of the config file to generate. + state, # A dictionary of state variables to use + logfile, # File handle where to write progress/logging output + particle_equil_events, # How many events will be run per particle to equilibrate the config. Useful if in setup you also need to equilibrate an intermediate configuration. +): + from subprocess import check_call + + #Here we work out how many unit cells to make the system out of for various packings + state = dict(state) + if 'InitState' not in state: + state['InitState'] = "FCC" + + unitcellN = { + "FCC":4, + "HCP":4, + "BCC":2, + "SC":1, + } + Ncells_unrounded = (state['N'] / unitcellN[state['InitState']]) ** (1.0 / 3.0) + Ncells = int(round(Ncells_unrounded)) + if abs(Ncells - Ncells_unrounded) > 0.1: + raise RuntimeError("Could not make "+str(state['N'])+" particles in an "+state['InitState']+" packing") + + # Here, for tethered systems, we do not simulate state points if + # its going to be boring and "ideal". I only have worked out the + # spacing expression for FCC, so all other crystals will just be + # run regardless + + if ("Rso" in state) and (state['Rso'] != float('inf')) and (state['InitState'] == "FCC"): + effrho = state['ndensity']*(state['Lambda']**3) + minR = max(0, (2**(2.5)*effrho)**(-1/3.0) - 0.5) + #phiT= state['ndensity'] * (4/3) * math.pi * minR**3 + #minRho = max(0, (2**(1/6.0)-(6*state['ndensity']*(4/3)*minR**3)**(1/3))**3) + if state['Rso'] <= minR: + raise pydynamo.SkipThisPoint() + + # This check is halting systems deep in the solid region, which should not be done! + # + # + ### Again, for tethered systems in FCC lattices we do not simulate + ### much beyond a multiple of the minimum tether radius. + ##if ("Rso" in state) and (state['ndensity'] >= 1.0) and (state['InitState'] == "FCC"): + ## minR = max(0, (2**(2.5)*state['ndensity'])**(-1/3.0) - 0.5) + ## if state['Rso'] >= 10*minR: + ## raise pydynamo.SkipThisPoint() + + + # Thermostat + options = '' + if 'kT' in state: + if state['kT'] != float('inf'): + options = options + ' -T '+repr(state['kT']) + else:#infite temperature is a special case, we set well energies to zero + options = options + ' -T 1.0' + + # Crystal lattice packing + packmode = { + 'FCC':0, + 'BCC':1, + 'SC': 2, + 'HCP':3, + } + options = options + ' --i1 '+str(packmode[state['InitState']]) +' -C '+str(Ncells) + + # Square well or hard sphere? + if state['Lambda'] != float('inf'): + if state['kT'] != float('inf'): + options = options + ' -m 1 --f1 '+repr(state['Lambda']) + else: # infinite temperature is a special case, we set well energies to zero + options = options + ' -m 1 --f1 '+repr(state['Lambda']) + " --f2 0.0" + else: + options = options + ' -m 0' + + # denisty + options = options + ' -d ' + repr(state['ndensity']) + + # Execution of dynamod + print('# dynamod'+options+' -o '+config, file=logfile) + check_call(('dynamod'+options+' -o '+config).split(), stdout=logfile, stderr=logfile) + + # Run of an equilibration step to blur the system state + if ('Rso' in state) and (state['Rso'] != float('inf')) and (state["InitState"] == "Liquid"): + print("\n", file=logfile) + print("################################", file=logfile) + print("# Liquifaction Run #", file=logfile) + print("################################\n", file=logfile, flush=True) + print("# dynarun --unwrapped "+config+" -o "+config+" -c "+str(state['N'] * particle_equil_events)+" --out-data-file data.liqequil.xml.bz2", file=logfile) + check_call(["dynarun", "--unwrapped", config, '-o', config, '-c', str(state['N'] * particle_equil_events), "--out-data-file", "data.liqequil.xml.bz2"], stdout=logfile, stderr=logfile) + + # Add the SO Cells global interaction (if needed) + if ('Rso' in state) and (state['Rso'] != float('inf')): + xml = pydynamo.ConfigFile(config) + XMLGlobals = xml.tree.find(".//Globals") + XMLSOCells = ET.SubElement(XMLGlobals, 'Global') + XMLSOCells.attrib['Name'] = "SOCells" #Name can be anything + XMLSOCells.attrib['Type'] = "SOCells" #This must be the right type of Global to load + XMLSOCellsRange = ET.SubElement(XMLSOCells, 'Range') + XMLSOCellsRange.attrib["Type"] = "All" + XMLSOCells.attrib['Diameter'] = str(2 * state['Rso']) + xml.save(config) + + +################################################################ +### DEFINE THE "STATE" VARIABLES TO BE SWEPT & RANGE +################################################################ +# This is the list of state variables and their ranges + + +statevars = [ + [ #Sweep + ("Lambda", [float('inf')]), + ("InitState", ["FCC"]), + ("N", list(map(lambda x: 4*x**3, [3,4,5,6,7]))), #15 + ('ndensity', list(set(map(lambda x : pydynamo.roundSF(x, 3), [0.01, 0.1, 0.5, 1.0, 1.3])))), + ("kT", [1.0]), + ], +] + +################################################################ +### CREATE A SIMULATION MANAGER +################################################################ +mgr = pydynamo.SimManager("HS_stats", #Which subdirectory to work in + statevars, #State variables + ["p", "CollisionMatrix"], # 'RadialDist' "VACF", # Output properties + restarts=1, #How many restarts (new initial configurations) should be done per state point + processes=None, #None is automatically use all processors +) + +################################################################ +### REORGANISE ANY EXISTING SIMULATIONS +################################################################ +#mgr.reorg_dirs() + +################################################################ +### RUN SOME SIMULATIONS +################################################################ +mgr.run(setup_worker=setup_worker, + particle_equil_events = 1000, # How many events per particle to equilibrate each sim for + particle_run_events = 1000, # How many events per particle to run IN TOTAL + particle_run_events_block_size=1000) # How big a block each run should be (for jacknife averaging). + +################################################################ +### GET THE DATA +################################################################ +# This creates a pandas dataframe with columns for the state variables +# AND any output values. It also generates pkl files, some for +# different properties. +df, state_data = mgr.fetch_data(1000) + +print(df) diff --git a/src/dynamo/dynamo/outputplugins/collMatrix.cpp b/src/dynamo/dynamo/outputplugins/collMatrix.cpp index 72f0dd1dd..ffa463622 100644 --- a/src/dynamo/dynamo/outputplugins/collMatrix.cpp +++ b/src/dynamo/dynamo/outputplugins/collMatrix.cpp @@ -32,6 +32,8 @@ void OPCollMatrix::initialise() { lastEvent.resize(Sim->N(), lastEventData(Sim->systemTime, EventKey(EventSourceKey(0, NOSOURCE), NONE))); + + _sysLastEventData = lastEventData(Sim->systemTime, EventKey(EventSourceKey(0, NOSOURCE), NONE)); // Reset the current capture state _currentCaptureState.clear(); @@ -174,7 +176,23 @@ void OPCollMatrix::eventUpdate(const Event &event, const NEventData &SDat) { newEvent(pData.particle1_.getParticleID(), pData.getType(), ck); newEvent(pData.particle2_.getParticleID(), pData.getType(), ck); } + + // Update the last event for the system as a whole + EventKey ek(ck, event._type); + if (_sysLastEventData.second.first.second != NOSOURCE) { + InterEventKey sysKey(ek, _sysLastEventData.second); + double dt = Sim->systemTime - _sysLastEventData.first; + //Perform an insert if the key doesn't exist, otherwise return the existing value + auto it = _sysInterEventMFTHistograms.insert(decltype(_sysInterEventMFTHistograms)::value_type( + sysKey, SysMFTData(Sim->lastRunMFT * 0.01 / Sim->N()))); + it.first->second.MFT.addVal(dt); + } + + _sysLastEventData.first = Sim->systemTime; + _sysLastEventData.second = ek; } + + void OPCollMatrix::newEvent(const size_t &part, const EEventType &etype, const EventSourceKey &ck) { if (lastEvent[part].second.first.second != NOSOURCE) { @@ -193,8 +211,26 @@ void OPCollMatrix::newEvent(const size_t &part, const EEventType &etype, void OPCollMatrix::output(magnet::xml::XmlStream &XML) { - XML << magnet::xml::tag("CollCounters") - << magnet::xml::tag("TransitionMatrix"); + XML << magnet::xml::tag("CollCounters"); + + XML << magnet::xml::tag("SystemMFT"); + + for (const auto &pair : _sysInterEventMFTHistograms) { + XML << magnet::xml::tag("MFT") << magnet::xml::attr("Event") + << pair.first.first.second << magnet::xml::attr("Name") + << getEventSourceName(pair.first.first.first, Sim) + << magnet::xml::attr("lastEvent") << pair.first.second.second + << magnet::xml::attr("lastName") + << getEventSourceName(pair.first.second.first, Sim); + + pair.second.MFT.outputHistogram(XML, 1.0 / Sim->units.unitTime()); + + XML << magnet::xml::endtag("MFT"); + } + + XML << magnet::xml::endtag("SystemMFT"); + + XML << magnet::xml::tag("TransitionMatrix"); std::map> totmap; @@ -274,8 +310,8 @@ void OPCollMatrix::output(magnet::xml::XmlStream &XML) { XML << magnet::xml::endtag("RijDotVij"); XML << magnet::xml::tag("RijDotDeltaPij"); - val.second.rijdotvij.outputHistogram(XML, 1.0 / Sim->units.unitLength() / - Sim->units.unitMomentum()); + val.second.rijdotdP.outputHistogram(XML, 1.0 / Sim->units.unitLength() / + Sim->units.unitMomentum()); XML << magnet::xml::endtag("RijDotDeltaPij"); XML << magnet::xml::tag("V2"); diff --git a/src/dynamo/dynamo/outputplugins/collMatrix.hpp b/src/dynamo/dynamo/outputplugins/collMatrix.hpp index c4e877806..3087b36f4 100644 --- a/src/dynamo/dynamo/outputplugins/collMatrix.hpp +++ b/src/dynamo/dynamo/outputplugins/collMatrix.hpp @@ -52,7 +52,7 @@ class OPCollMatrix : public OutputPlugin { unsigned long totalCount; // EventKey is a pair of EventSourceKey and EEventType - // It describes the type of event and its source + // Combines the EventSourceKey (type i.e. INTERACTION, and ID) with the EEventType (CORE, WELL, WALL, VIRTUAL, etc) //! A key for two events typedef std::pair InterEventKey; @@ -110,8 +110,21 @@ class OPCollMatrix : public OutputPlugin { std::map initialCounter; + //! The time and event key for the last event for each particle typedef std::pair lastEventData; + //! Keeps track of the last event for each particle std::vector lastEvent; + + //! Keeps track of the last event for the system as a whole + lastEventData _sysLastEventData; + + struct SysMFTData { + SysMFTData(double binWidth) + : MFT(binWidth) {} + magnet::math::Histogram<> MFT; + }; + //! Keeps track of the MFT histograms for the system as a whole + std::map _sysInterEventMFTHistograms; }; } // namespace dynamo diff --git a/src/dynamo/dynamo/outputplugins/eventtypetracking.hpp b/src/dynamo/dynamo/outputplugins/eventtypetracking.hpp index 11f1512de..7e8a38638 100644 --- a/src/dynamo/dynamo/outputplugins/eventtypetracking.hpp +++ b/src/dynamo/dynamo/outputplugins/eventtypetracking.hpp @@ -27,10 +27,10 @@ class System; namespace EventTypeTracking { -//! Keeps the ID and type of the event source +//! Keeps the type of the event source (GLOBAL, LOCAL, INTERACTION, SYSTEM) and the ID typedef std::pair EventSourceKey; -//! Event source And Type +//! Combines the EventSourceKey (type i.e. INTERACTION, and ID) with the EEventType (CORE, WELL, WALL, VIRTUAL, etc) typedef std::pair EventKey; std::string getEventSourceName(const EventSourceKey &, From 59f4bb54a0a3f8c2f2e71c58597d1c2e75922d2d Mon Sep 17 00:00:00 2001 From: Marcus Bannerman Date: Fri, 4 Sep 2026 14:16:53 +0100 Subject: [PATCH 18/19] Parsing the SysMFT --- src/pydynamo/output_properties.py | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/src/pydynamo/output_properties.py b/src/pydynamo/output_properties.py index 53ac838ce..78d686a1e 100644 --- a/src/pydynamo/output_properties.py +++ b/src/pydynamo/output_properties.py @@ -346,7 +346,15 @@ def result(self, state, outputfile: OutputFile, configfilename, counter, manager captures = int(tag.attrib["captures"]) # Square velocity retval["V2"][(iName, eType, captures)] += Histogram.load(tag.find("./V2/Histogram")) - + + for tag in outputfile.tree.findall(".//CollCounters/SystemMFT/MFT"): + iName1 = tag.attrib["Name"] + eType1 = tag.attrib["Event"] + iName2 = tag.attrib["lastName"] + eType2 = tag.attrib["lastEvent"] + # MFT histogram + retval["SysMFT"][(iName1, eType1, iName2, eType2)] += Histogram.load(tag.find("./Histogram")) + return retval def init(self): @@ -354,6 +362,7 @@ def init(self): "CollisionMatrix": KeyedWeightedKeyedArray(), "CaptureStateHistogram": Histogram(), "V2": KeyedHistogram(), + "SysMFT": KeyedHistogram(), } OutputFile.output_props["CollisionMatrix"] = CollisionMatrixOutputProperty() From 3eeeffade2271f04fece4736852e2fa4076f2b23 Mon Sep 17 00:00:00 2001 From: Marcus Bannerman Date: Fri, 4 Sep 2026 14:34:08 +0100 Subject: [PATCH 19/19] Added ignoring virtual events --- .../dynamo/outputplugins/collMatrix.cpp | 25 +++++++++++-------- 1 file changed, 15 insertions(+), 10 deletions(-) diff --git a/src/dynamo/dynamo/outputplugins/collMatrix.cpp b/src/dynamo/dynamo/outputplugins/collMatrix.cpp index ffa463622..d02f2bce4 100644 --- a/src/dynamo/dynamo/outputplugins/collMatrix.cpp +++ b/src/dynamo/dynamo/outputplugins/collMatrix.cpp @@ -179,17 +179,22 @@ void OPCollMatrix::eventUpdate(const Event &event, const NEventData &SDat) { // Update the last event for the system as a whole EventKey ek(ck, event._type); - if (_sysLastEventData.second.first.second != NOSOURCE) { - InterEventKey sysKey(ek, _sysLastEventData.second); - double dt = Sim->systemTime - _sysLastEventData.first; - //Perform an insert if the key doesn't exist, otherwise return the existing value - auto it = _sysInterEventMFTHistograms.insert(decltype(_sysInterEventMFTHistograms)::value_type( - sysKey, SysMFTData(Sim->lastRunMFT * 0.01 / Sim->N()))); - it.first->second.MFT.addVal(dt); - } - _sysLastEventData.first = Sim->systemTime; - _sysLastEventData.second = ek; + // We ignore virtual events for the system as a whole, since they don't represent dynamics + if (event._type != VIRTUAL) { + // Ignore the first event, since we don't have a previous event to compare to + if (_sysLastEventData.second.first.second != NOSOURCE) { + InterEventKey sysKey(ek, _sysLastEventData.second); + double dt = Sim->systemTime - _sysLastEventData.first; + //Perform an insert if the key doesn't exist, otherwise return the existing value + auto it = _sysInterEventMFTHistograms.insert(decltype(_sysInterEventMFTHistograms)::value_type( + sysKey, SysMFTData(Sim->lastRunMFT * 0.01 / Sim->N()))); + it.first->second.MFT.addVal(dt); + } + + _sysLastEventData.first = Sim->systemTime; + _sysLastEventData.second = ek; + } }