diff --git a/.gitignore b/.gitignore index 9f98cc17f..a58805ca4 100644 --- a/.gitignore +++ b/.gitignore @@ -12,4 +12,9 @@ __pycache__ Testing/Temporary/CTestCostData.txt .eggs wheelhouse -vcpkg_installed \ No newline at end of file +vcpkg_installed + +result +build_* +result +result-* \ No newline at end of file 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/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 1ff59a038..90d802990 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 ? false, + }: 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 + bzip2 + boost eigen - + ] ++ (lib.optionals visualiser [ # Visualiser libGL gtkmm3.dev @@ -47,5 +51,5 @@ python3.pkgs.buildPythonPackage rec { cairomm.dev libpng mesa - ]; + ]); } 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/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..f590192ff --- /dev/null +++ b/scripts/plotter.ipynb @@ -0,0 +1,108 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "id": "bbb1036b", + "metadata": {}, + "outputs": [], + "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", + "# 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": null, + "id": "3c5f65e0", + "metadata": {}, + "outputs": [], + "source": [ + "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", + "\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", + " 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.yscale(\"log\")\n", + "plt.legend()\n", + "plt.grid()\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fc6e9d0c", + "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/dynamo/dynamo/outputplugins/collMatrix.cpp b/src/dynamo/dynamo/outputplugins/collMatrix.cpp index f36497c4b..d02f2bce4 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(); @@ -76,8 +78,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); @@ -160,7 +176,28 @@ 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); + + // 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; + } } + + void OPCollMatrix::newEvent(const size_t &part, const EEventType &etype, const EventSourceKey &ck) { if (lastEvent[part].second.first.second != NOSOURCE) { @@ -179,8 +216,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; @@ -260,8 +315,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"); @@ -271,8 +326,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..3087b36f4 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 + // EventKey is a pair of EventSourceKey and EEventType + // Combines the EventSourceKey (type i.e. INTERACTION, and ID) with the EEventType (CORE, WELL, WALL, VIRTUAL, etc) //! 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,17 +87,44 @@ 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; 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 &, diff --git a/src/pydynamo/__init__.py b/src/pydynamo/__init__.py index b7c9a8db9..ba7429e68 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): @@ -185,12 +185,21 @@ def perdir(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 @@ -313,7 +322,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 @@ -518,13 +527,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 +564,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 +573,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 diff --git a/src/pydynamo/output_properties.py b/src/pydynamo/output_properties.py index 7d10d2f63..78d686a1e 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, KeyedKeyedArray, + WeightedType, Histogram, KeyedHistogram) # A XMLFile/ElementTree but specialised for DynamO output files @@ -48,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 @@ -65,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) @@ -104,7 +106,23 @@ 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)) + 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 + [], # 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 = [] @@ -117,9 +135,7 @@ 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): def __init__(self): @@ -134,14 +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 @@ -181,7 +199,8 @@ 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): def __init__(self): @@ -195,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) @@ -208,7 +227,8 @@ 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): ''' @@ -217,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 @@ -230,26 +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) + return {"FCCOrder":WeightedType()} +OutputFile.output_props["OrderParameter"] = OrderParameterProperty() class ChungLuConfigurationModel(OutputProperty): ''' @@ -261,7 +270,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 +284,85 @@ 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["ChungLu"]["bond_order_count"] = WeightedType(bond_order_counter, 1) + retval["ChungLu"]["order_count"] = WeightedType(order_count, 1) + return retval def init(self): - return WeightedType(KeyedArray(), 0) + 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() + + # 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"] + 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()) + + # 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")) + + 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): + return { + "CollisionMatrix": KeyedWeightedKeyedArray(), + "CaptureStateHistogram": Histogram(), + "V2": KeyedHistogram(), + "SysMFT": KeyedHistogram(), + } -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 6f6f25335..8527b84ed 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) @@ -96,4 +96,57 @@ 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"]["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(float, {"1":2}) + + c = a + b + 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"]["1"] == pytest.approx(-1) + assert c["b"]["2"] == pytest.approx(3) + assert c["c"]["3"] == pytest.approx(4) + + c = a * b + assert c["a"]["1"] == pytest.approx(2) + assert "b" not in c + 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 86ec02b92..5af2ad81b 100644 --- a/src/pydynamo/weighted_types.py +++ b/src/pydynamo/weighted_types.py @@ -10,18 +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. """ - - store = defaultdict(float) - - def __init__(self, values = None): - self.store = defaultdict(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") @@ -116,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): @@ -137,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) @@ -152,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) @@ -166,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: @@ -188,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) @@ -226,7 +225,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) @@ -236,4 +235,119 @@ 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, weight = 0.0): + super().__init__(value= KeyedArray() if (values == None) else values, weight=weight) + +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): + self.binwidth = binwidth + self.data = WeightedKeyedArray() + + def __repr__(self): + return f"Histogram(binwidth={self.binwidth}, 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) + new_histogram.data = self.data + other.data + return new_histogram + + def avg(self): + # 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): + """ + 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"]) + + 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())) + 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) + + return retval + +class KeyedHistogram(KeyedArray): + def __init__(self): + super().__init__(type=Histogram) + from collections import defaultdict