Speed between pyobjcryst-Debian and FOX-GUI-Windows

[hubsilon]

Dear all,

I hope this email finds you well.

I've been running the cimetadine example through pyobjcryst on debian via python scrypting (https://www.diffpy.org/pyobjcryst/examples/structure-solution-powder-cimetidine.html) and FOX on windows via graphical interface (starting from the xmlgz obtained from pyobjcryst in order to insure the same initial parameters)

I launched 3 runs of 100000 iterations in both cases but I observed some difference in calculation time dedicated to the global optimisation: 140 secondes for pyobjcryst and 55 secondes for FOX.

I'm thus wondering:

1. if I missed something with pyobjcryst (memory usage, etc...)?
2. if it is coming from python vs C++? ...but my guess is that python is only acting as an API while the real calcution are runned with C++)
3. if there is a way to accelerate it? ...by reducing time not parallelisation

thank you in advance for your advices

best regards

[vincefn]

One important difference between Fox and pyobjcryst is that Fox uses single precision, and pyobjcryst double. Though with modern processors I don't expect a huge difference (more like 1.5x ?).

Besides this I do not know - you certainly want to disable all plots in python (if you use a notebook), as you may lose a little time updating the plots (in Fox this is multithreaded).

I'd also double-check the extent (max sin theta / lambda) but if you load the xmlgz it should be the same.

Otherwise there should not be any significant difference - indeed everything is happening in c++, except for plotting if any.

[hubsilon]

Thank you for your quick reply.

I disable all display, only print a few parameters to insure it's running correctly.

I added a profiling on the part dedicated to Monte Carlo Optimisation: mc.MultiRunOptimize(nb_run=3, nb_step=1000000)

and obtained this result:
_LSQ option: Every 150000 trials, and at the end of each run
=== Profiling of Monte Carlo ===
=== End of profiling ===

     123 function calls in 138.519 seconds

Ordered by: cumulative time
ncalls tottime percall cumtime percall filename:lineno(function)
1 138.517 138.517 138.519 138.519 /tmp/sshfs-persistent_area/miniconda/lib/python3.12/site-packages/pyobjcryst/globaloptim.py:51(MultiRunOptimize)
40 0.001 0.000 0.001 0.000 /tmp/sshfs-persistent_area/miniconda/lib/python3.12/site-packages/pyobjcryst/globaloptim.py:122(UpdateDisplay)
40 0.000 0.000 0.000 0.000 /tmp/sshfs-persistent_area/miniconda/lib/python3.12/site-packages/pyobjcryst/crystal.py:73(UpdateDisplay)
40 0.000 0.000 0.000 0.000 /tmp/sshfs-persistent_area/sdpd/50_fox/cimetadine/diffpy_example.py:114()
1 0.000 0.000 0.000 0.000 /tmp/sshfs-persistent_area/miniconda/lib/python3.12/site-packages/pyobjcryst/globaloptim.py:86(_fix_parameters_for_global_optim)
1 0.000 0.000 0.000 0.000 {method 'disable' of 'lsprof.Profiler' objects}

With FOX on windows, it's more complicated for me, I used a stopwatch ;)

In order to give you a better idea, I attached below my pyobjcryst code which is a copy/paste of the cimetadine example running on a debian virtual machine with 16CPU and 64GB of memory:

import os
import pyobjcryst
import numpy as np
from pyobjcryst.crystal import *
from pyobjcryst.powderpattern import *
from pyobjcryst.indexing import *
from pyobjcryst.molecule import *
from pyobjcryst.globaloptim import MonteCarlo
from pyobjcryst.io import xml_cryst_file_save_global
import cProfile, pstats, io

p = PowderPattern()

data_dir = "./data"
dat_path = os.path.join(data_dir, "cime.dat")
if not os.path.exists(dat_path):
    os.makedirs(data_dir, exist_ok=True)
    os.system("curl -o {} https://raw.githubusercontent.com/vincefn/objcryst/master/Fox/example/tutorial-cimetidine/cime.dat".format(dat_path))

p.ImportPowderPatternFullprof(dat_path)
p.SetWavelength(1.52904)


# Index
pl = p.FindPeaks(1.5, -1, 1000)
if len(pl) > 20:
    pl.resize(20)  # Only keep 20 peaks
for peak in pl:
    print(peak)

ex = quick_index(pl)

print("Solutions:")
for s in ex.GetSolutions():
    print(s)
      
    
uc = ex.GetSolutions()[0][0].DirectUnitCell()
c = pyobjcryst.crystal.Crystal(uc[0], uc[1], uc[2], uc[3], uc[4], uc[5], "P1")
pdiff = p.AddPowderPatternDiffraction(c)

# Plot with indexing in new figure
#p.plot(diff=False,fig=None,hkl=True)


p.SetMaxSinThetaOvLambda(0.3)
p.quick_fit_profile(auto_background=True,plot=False, init_profile=True,verbose=True)
p.quick_fit_profile(plot=False, init_profile=False, asym=True, displ_transl=True, verbose=False)

# Plot in new figure
#p.plot(diff=True, fig=None, hkl=True)
print("Fit result: Rw=%6.2f%% Chi2=%10.2f  GoF=%8.2f  LLK=%10.3f" %
      (p.rw * 100, p.chi2, p.chi2/p.GetNbPointUsed(), p.llk))
      
      
      
      
p.SetMaxSinThetaOvLambda(0.2)  # Important for stability of profile fit. And faster !
spgex = SpaceGroupExplorer(pdiff)

# NB:verbose C++ output does not appear in a notebook
spgex.RunAll(keep_best=True, update_display=False, fitprofile_p1=False)

for sol in spgex.GetScores():
    #if sol.nGoF > 4 * spgex.GetScores()[0].nGoF:
    if sol.GoF <= 2 * spgex.GetScores()[0].GoF:
        print(sol)

print("Chosen spacegroup (smallest nGoF): ", c.GetSpaceGroup())

# Updated plot with optimal spacegroup
#p.plot(diff=True, fig=None, hkl=True, reset=True)



fhz_path = os.path.join(data_dir, "cime.fhz")
if not os.path.exists(fhz_path):
    os.system("curl -o {} https://raw.githubusercontent.com/vincefn/objcryst/master/Fox/example/tutorial-cimetidine/cime.fhz".format(fhz_path))

# This will automatically add the Molecule to the Crystal object
m = ImportFenskeHallZMatrix(c, fhz_path)
print("Crystal Formula:", c.GetFormula())
print()

# List the options to see which are available
for i in range(c.GetNbOption()):
    o = c.GetOption(i)
    print("Option #%d: %s = %s"% (i, o.GetName(), o.GetChoiceName(o.GetChoice())))
    for j in range (o.GetNbChoice()):
        print("   Choice %d: %s" %(j, o.GetChoiceName(j)))


c.GetOption(1).SetChoice(0)

mc = MonteCarlo()
mc.AddRefinableObj(c)
mc.AddRefinableObj(p)

pdiff.SetExtractionMode(False)
p.FitScaleFactorForRw() # for a better initial display
#p.plot(fig=None,diff=False,hkl=True)

mc.GetOption("Automatic Least Squares Refinement").SetChoice(2)
print("LSQ option: ", mc.GetOption("Automatic Least Squares Refinement").GetChoiceName(2))



# --- Profiling ---
print("\n=== Profiling of Monte Carlo ===\n")
pr = cProfile.Profile()
pr.enable()

# Deactivate refreshing of display
p.UpdateDisplay = lambda *args, **kwargs: None
p.plot = lambda *args, **kwargs: None

# launch MC: 3 runs of 1e6 iter
mc.MultiRunOptimize(nb_run=3, nb_step=1000000)

pr.disable()
print("\n=== End of profiling ===\n")

# save profiling report
s = io.StringIO()
sortby = pstats.SortKey.CUMULATIVE  
ps = pstats.Stats(pr, stream=s).sort_stats(sortby)
ps.print_stats(30)  
print(s.getvalue())

print("Final LLK: %.2f" % mc.GetLogLikelihood())


# Save result so it can be opened by Fox
xml_cryst_file_save_global('result.xmlgz')
# Also export to the CIF format
c.CIFOutput("result.cif")

I'd be keen to decrease the calculation time down to those with FOX or even better would be great.

[vincefn]

I tested on my macbook with a similar config - taking the example cimetidine-powder.xml, activating least squares every 150k, and running 3 times for 1M trials.

I got almost identical timings from Fox (no graphs) and running it in a notebook using pyobjcryst - both around 55s.

Here's the simplified code (I also read the objects from the xml file):

import os
import pyobjcryst
import numpy as np
from pyobjcryst.crystal import *
from pyobjcryst.powderpattern import *
from pyobjcryst.indexing import *
from pyobjcryst.molecule import *
from pyobjcryst.globaloptim import MonteCarlo
from pyobjcryst.io import xml_cryst_file_save_global, xml_cryst_file_load_all_object
from pyobjcryst._pyobjcryst import gOptimizationObjRegistry
import cProfile, pstats, io

c,p =xml_cryst_file_load_all_object("/Users/favre/Desktop/CimeTest.xmlgz")
mc = gOptimizationObjRegistry.GetObj(0)

mc.GetOption("Automatic Least Squares Refinement").SetChoice(2)
print("LSQ option: ", mc.GetOption("Automatic Least Squares Refinement").GetChoiceName(2))



# --- Profiling ---
print("\n=== Profiling of Monte Carlo ===\n")
pr = cProfile.Profile()
pr.enable()

# Deactivate refreshing of display
p.UpdateDisplay = lambda *args, **kwargs: None
p.plot = lambda *args, **kwargs: None

# launch MC: 3 runs of 1e6 iter
mc.MultiRunOptimize(nb_run=3, nb_step=1000000)

pr.disable()
print("\n=== End of profiling ===\n")

# save profiling report
s = io.StringIO()
sortby = pstats.SortKey.CUMULATIVE  
ps = pstats.Stats(pr, stream=s).sort_stats(sortby)
ps.print_stats(30)  
print(s.getvalue())

print("Final LLK: %.2f" % mc.GetLogLikelihood())


# Save result so it can be opened by Fox
xml_cryst_file_save_global('result.xmlgz')

You are running under a debian virtual machine - so this is not the same computer/processor than under windows, so hard to make a comparison. Many-cores machines can often have slower individual core speed, whereas desktop machines can allow single-core runs to be boosted.

Same test under linux/ubuntu 24.04 (AMD EPYC 9354 32-Core @3.8 GHz) takes 72s in Fox and 94s in pyobjcryst. I expect the difference here is due to the double precision in pyobjcryst I mentioned earlier.