regmaxs/regmaxsn/core/transOnce.py

from SWCTransforms import SWCTranslate, ArgGenIterator, objFun
import multiprocessing as mp
import numpy as np
import json
import sys
from itertools import product
from regmaxsn.core.transforms import compose_matrix

debugging = False
# debugging = True

assert len(sys.argv) == 2, 'Only one argument, the path of the swcfile expected, ' + str(len(sys.argv)) + 'found'

parFile = sys.argv[1]
with open(parFile, 'r') as fle:
    pars = json.load(fle)

refSWC, SWC2Align, outFiles, gridSizes, bounds, minRes, nCPU = pars

SWCDatas = [SWCTranslate(refSWC, SWC2Align, x) for x in gridSizes]
pool = mp.Pool(processes=nCPU)
bestSol = [0, 0, 0]

for gridInd, gridSize in enumerate(gridSizes):


    bounds = (np.array(bounds).T - np.array(bestSol)).T
    boundsRoundedUp = np.sign(bounds) * np.ceil(np.abs(bounds) / gridSize) * gridSize
    possiblePts1D = [(bestSol[ind] + np.arange(x[0], x[1] + gridSize, gridSize)).tolist()
                     for ind, x in enumerate(boundsRoundedUp)]

    possiblePts3D = list(product(*possiblePts1D))
    if debugging:
        print('Gridsize:' + str(gridSize))
        print(bounds)
        print(map(len, possiblePts1D))
        print([bestSol[ind] + x for ind, x in enumerate(boundsRoundedUp)])

    argGen = ArgGenIterator(possiblePts3D, SWCDatas[gridInd])
    funcVals = pool.map_async(objFun, argGen).get(1800)
    minimum = min(funcVals)
    minimzers = [y for x, y in enumerate(possiblePts3D) if funcVals[x] == minimum]

    if not gridInd:
        distFrom0 = np.linalg.norm(minimzers, axis=1)
        bestSol = minimzers[np.argmin(distFrom0)]
    else:

        prevVals = [objFun((x, SWCDatas[gridInd - 1])) for x in minimzers]
        bestSol = minimzers[np.argmin(prevVals)]
    bounds = map(lambda x: [x - gridSize, x + gridSize], bestSol)
    if debugging:
        bestVal = objFun((bestSol, SWCDatas[gridInd]))
        print(bestSol, bestVal)

if minRes < gridSizes[-1]:

    bounds = (np.array(bounds).T - np.array(bestSol)).T
    boundsRoundedUp = np.sign(bounds) * np.ceil(np.abs(bounds) / minRes) * minRes
    possiblePts1D = [(bestSol[ind] + np.arange(x[0], x[1] + minRes, minRes)).tolist()
                     for ind, x in enumerate(boundsRoundedUp)]
    possiblePts3D = list(product(*possiblePts1D))

    if debugging:
        print('StepSize:' + str(minRes))
        print(bounds)
        print(map(len, possiblePts1D))
        print([bestSol[ind] + x for ind, x in enumerate(boundsRoundedUp)])

    argGen = ArgGenIterator(possiblePts3D, SWCDatas[-1])
    funcVals = pool.map_async(objFun, argGen).get(1800)

    minimum = min(funcVals)

    minimzers = [y for x, y in enumerate(possiblePts3D) if funcVals[x] == minimum]
    prevVals = [objFun((x, SWCDatas[-2])) for x in minimzers]
    bestSol = minimzers[np.argmin(prevVals)]

bestVal = objFun((bestSol, SWCDatas[-1]))
nochange = objFun(([0, 0, 0], SWCDatas[-1]))
if debugging:
    bestVals = [objFun((bestSol, x)) for x in SWCDatas]
    print(bestSol, bestVals)

done = False

# all values are worse than doing nothing
if bestVal > nochange:

    done = True
    bestSol = [0, 0, 0]
    bestVal = nochange

# best solution and no change are equally worse
elif bestVal == nochange:

    # the solution step is very close to zero
    if np.abs(bestSol).max() <= min(minRes, gridSizes[-1]):

        done = True
        bestSol = [0, 0, 0]
        bestVal = nochange



SWCDatas[-1].writeSolution(outFiles[0], bestSol)
matrix = compose_matrix(translate=bestSol).tolist()
with open(outFiles[1], 'w') as fle:
    json.dump({'type': 'XYZ Translations in um','bestSol': bestSol,
               'transMat': matrix, 'done': done, 'bestVal': bestVal}, fle)