############################################################################################################################
#
# Data load setup
#
############################################################################################################################

# img out sample calls
#
# gen blended sequence for 2 inputs:
#./manta ../mantaflowgit/scenes/pcgTestOfLoad.py doOnlyLoad 1 doSingleLoad 1 dataSet 201 showUi 1 blendIters 1 blendSingleAlpha 0.5 regen2dSlices 1 dispWidth 500 dispHeight 500 camposz 0.8 singleSel   0 0 0   0 1 0
#
# gen blended sequence for 3 inputs:
#./manta ../mantaflowgit/scenes/pcgTestOfLoad.py doOnlyLoad 1 doSingleLoad 1 dataSet 202 showUi 1 blendIters 1 regen2dSlices 1 dispWidth 500 dispHeight 500 camposz 0.8 singleSel3   0 0 0   0 1 0   1 0 0      blendSingleAlpha 0.5   blendThirdAlpha 0.5
#
# run single grid data
#./manta ../mantaflowgit/scenes/pcgTestOfLoad.py dataSet 202  dosingleselect 1  singleSel 0 1 0  0 2 0
#

import sys
from manta import *
from math  import *

setDebugLevel(1); # set to 2 for kernel msgs, 3 for ranges


dataSet       = 2
paramSet      = 5
showUi        = True
multiStep     = 3
doLsProject   = False
doFinalProject= True
doStoreResult = False
writePngs     = False

# load only modes
doOnlyLoad     = False
doSingleSelect = False
singleSel      = [0,0,0, 1,1,0,  -1,-1,-1]; 
doLoadThird    = False
doConstOverride = 0; # for testing only, paper ex2, turn of SDF!

# configure blending
blendIters       = -9;
blendSingleAlpha = -1;
blendThirdAlpha = -1;

# generate 2d slices from 3d version
regen2dSlices = False
# window settings, for img output
dispWidth  = 800
dispHeight = 600
camPosZ    = 1.2

# data input 
dim = 2
res = 10
resLoad        = 12
isLevelset     = True
frameIterStart = 0
frameIterEnd   = 12
frameIterMult  = 10
frameIterOffset= 10
orderTime  = 1
orderSpace = 1
projSizeThreshold = 959
doLoadGrid = 0
thirdDimFactor = 1.0

# modify filename for some runs...
dataSetSuffix = ""

# fit data into grid with 10% border, used in OF solve for reset
autoBorder  = 0.1
loadOffset  = Vec3(0,0,0);

gridLoadX  = -1;  gridLoadY = -1;  gridLoadZ = -1;
gridOtherX = -1; gridOtherY = -1; gridOtherZ = -1;
gridThirdX = -1; gridThirdY = -1; gridThirdZ = -1;


# check for parameters (ignore case), note: bool given as 0/1
# eg call with: "./manta scene.py   dataSet 2 showUi 1 doOnlyLoad 0 "
skipNext = 0;
for iter in range(1, len(sys.argv)):
	if skipNext>0:
		skipNext = skipNext-1; continue; 
	if(sys.argv[iter].lower() == "showui"):
		showUi = (int(sys.argv[iter+1])!=0);
		print( "showUi set to %d" % (showUi) ); skipNext = 1;
	elif(sys.argv[iter].lower() == "dataset"):
		dataSet = int(sys.argv[iter+1])
		print( "dataSet set to %d" % (dataSet) ); skipNext = 1;
	elif(sys.argv[iter].lower() == "paramset"):
		paramSet = int(sys.argv[iter+1])
		print( "paramSet set to %d" % (paramSet) ); skipNext = 1;
	elif(sys.argv[iter].lower() == "dostoreresult"):
		doStoreResult = (int(sys.argv[iter+1])!=0);
		print( "doStoreResult set to %d" % (doStoreResult) ); skipNext = 1;
	elif(sys.argv[iter].lower() == "doonlyload"):
		doOnlyLoad = (int(sys.argv[iter+1])!=0);
		print( "doOnlyLoad set to %d" % (doOnlyLoad) ); skipNext = 1;
	elif(sys.argv[iter].lower() == "dosingleload") or (sys.argv[iter].lower() == "dosingleselect"):
		doSingleSelect = (int(sys.argv[iter+1])!=0);
		print( "doSingleSelect set to %d" % (doSingleSelect) ); skipNext = 1;
	elif(sys.argv[iter].lower() == "singlesel"):
		for i in range(0, 6):
			singleSel[i] = int(sys.argv[iter+1+i]);
		print( "singleSel set to %d,%d,%d and %d,%d,%d" % (singleSel[0],singleSel[1],singleSel[2],singleSel[3],singleSel[4],singleSel[5]) ); 
		skipNext = 6;
	elif(sys.argv[iter].lower() == "singlesel3"): # double interpol, 9 values
		for i in range(0, 9):
			singleSel[i] = int(sys.argv[iter+1+i]);
		print( "singleSel3 set to " + str(singleSel) ); 
		skipNext = 9;
	elif(sys.argv[iter].lower() == "blenditers"):
		blendIters = int(sys.argv[iter+1])
		print( "blendIters set to %d" % (blendIters) ); skipNext = 1;
	elif(sys.argv[iter].lower() == "blendsinglealpha"):
		blendSingleAlpha = float(sys.argv[iter+1])
		print( "blendSingleAlpha set to %f" % (blendSingleAlpha) ); skipNext = 1;
	elif(sys.argv[iter].lower() == "blendthirdalpha"):
		blendThirdAlpha = float(sys.argv[iter+1])
		print( "blendThirdAlpha set to %f" % (blendThirdAlpha) ); skipNext = 1;
	elif(sys.argv[iter].lower() == "regen2dslices"):
		regen2dSlices = (int(sys.argv[iter+1])!=0);
		print( "regen2dSlices set to %d" % (regen2dSlices) ); skipNext = 1;
	elif(sys.argv[iter].lower() == "dispwidth"):
		dispWidth = int(sys.argv[iter+1])
		print( "dispWidth set to %d" % (dispWidth) ); skipNext = 1;
	elif(sys.argv[iter].lower() == "dispheight"):
		dispHeight = int(sys.argv[iter+1])
		print( "dispHeight set to %d" % (dispHeight) ); skipNext = 1;
	elif(sys.argv[iter].lower() == "camposz"):
		camPosZ = float(sys.argv[iter+1])
		print( "camPosZ set to %d" % (camPosZ) ); skipNext = 1;
	else:
		print( "Error: unknown parameter '" + sys.argv[iter] +"' " );
		exit(1);



############################################################################################################################

# liquids

# 1-3 , 64^2 in pcgData64
if dataSet==1:
	datPath = "pcgData64/phi"
	frameIterEnd = 8
	res = 64
	dat0 = 3
	dat1 = 1  
	resLoad = 64
	#res = 100
	#frameIterOffset= 60; # more interesting start frame
	#dat0 = 1; dat1 = 2; # NT_DEBUG

# 4-6 , 140^2 in pcgData140
elif dataSet==2:
	datPath = "pcgData140/phi"
	res = 140
	resLoad = 140
	res = 168
	frameIterStart = 4 # debug , tough start frame

	# variants
	#dat0 = 4; dat1 = 5; # rel. easy
	dat0 = 4; dat1 = 6;  # tough!
	#dat0 = 6;	dat1 = 4;  # tougher , reverse
	#dat0 = 5; dat1 = 4; # NT_DEBUG new


# border (autoBord,autoBnd) reduce test data
elif dataSet==3:
	autoBorder  = 0.0
	datPath = "pcgDataBordred/phibv"
	res = 120
	dat0 = 1; 
	dat1 = 3; 
	resLoad = 100
	frameIterStart = 2 
	frameIterMult  = 1
	frameIterOffset= 1

# smoke

# 2d plume 1-4 (note 1&2 are too low...?)
elif dataSet==5:
	datPath = "pcgDataSmoke64/dens"
	#frameIterOffset = 5;
	frameIterOffset = 30;
	frameIterMult   = 5;
	frameIterEnd    = 8;
	res = 64;
	isLevelset = False;
	dat0 = 3;
	dat1 = 4;
	resLoad = 64;
	res = 80;
	frameIterEnd    = 18; frameIterOffset = 5 # longer , all frames
	#frameIterOffset = 60 # start later

# 2d paper example
elif dataSet==6:
	datPath = "pcgData2dPaper/phi"
	res = 140
	dat0 = 1; 
	dat1 = 2;
	resLoad = 128
	res = 280;
	frameIterStart = 2 

# 3d

# 7-9 , 64^3 in pcgData64 (3d versions of 1-3)
elif dataSet==10:
	datPath = "pcgData64threed/phi"
	dim = 3; res = 64
	dat0 = 7
	dat1 = 9
	resLoad = 64
	res = 16; # NT_DBEUG 4d
	# shorter
	frameIterStart = 1; frameIterEnd   = 7
	frameIterMult  = 20; frameIterOffset= 10
	#res = 124;	frameIterMult   = 5;	frameIterEnd    = 80; 


# 3d plume 5,6
elif dataSet==15:
	datPath = "pcgDataSmoke64threed/dens"
	frameIterOffset = 50
	frameIterMult   = 5
	dim = 3; res = 64
	isLevelset = False
	dat0 = 5
	dat1 = 6
	resLoad = 64
	# full run
	frameIterOffset = 5; frameIterEnd = 18

# 3d volumes for dat 1-3 , 64^2 in pcgData64conv3d , frame #0 is interval 1, #1 is interval 2
elif dataSet==101:
	datPath = "pcgData64conv3d/phito"
	dim = 3;
	res = 64;
	dat0 = 1;
	dat1 = 3;
	dat0 = 3; dat1 = 1; # reverse
	resLoad = 64;
	res = 76;
	frameIterOffset= 0;
	frameIterMult  = 1;
	frameIterEnd   = 2;
	#resLoad = 64; res = 40; # smaller
	frameIterOffset= 1; frameIterEnd   = 1;

elif dataSet==102:
	# special test, run a case from 202 single, for loop
	dat0 = 1; dat1 = 2; # fixed for grids!
	datPath = "pcgDataGrid01/phito02" 
	doLoadGrid = 1;
	dim        = 3;
	resLoad    = 80;
	res        = 96;
	frameIterOffset= 0; frameIterMult  = 1;
	frameIterStart = 0; frameIterEnd   = 1; 
	gridLoadX,  gridLoadY,  gridLoadZ  = (1,5,0)
	gridLoadX,  gridLoadY,  gridLoadZ  = (1,1,0) # tougher!
	gridOtherX, gridOtherY, gridOtherZ = (3,5,0)
	res        = 44; # small! 
	res        = 64; # ca. 10s

elif dataSet==103:
	# special test, run a case from 202 single, for loop
	dat0 = 1; dat1 = 2; # fixed for grids!
	datPath = "pcgDataFlip2d/phifl2d" 
	doLoadGrid = 0;
	dim        = 3;
	resLoad    = 100;
	res        = 120;
	frameIterOffset= 0; frameIterMult  = 1;
	frameIterStart = 0; frameIterEnd   = 1; 
	res        = 64; # small! 

elif dataSet==104:
	#NT_DEBUG print("ERR fix time below"); exit(1);
	# two d case with longer time
	dat0 = 1; dat1 = 2; 
	datPath = "pcgDataFlip2d/phifl2dLonger" 
	doLoadGrid = 0;
	dim        = 3;
	resLoad    = 50;
	res        = 60;
	frameIterOffset= 0; frameIterMult  = 1;
	frameIterStart = 0; frameIterEnd   = 1; 
	#res        = 34; # small! 
	thirdDimFactor = 1.5


# 3d volumes for 2d flip sims, grid 10x10;   to02: 80 , to03: 46  (pcgTestDatagenArrFlip)
elif dataSet==200:
	dat0 = 1; dat1 = 2; # fixed for grids!
	datPath = "pcgDataGrid01/phito03" # note, suffix added!
	doLoadGrid = 1;
	dim        = 3;
	resLoad    = 46;
	res        = 56;
	#res        = 32; # smaller NT DEBUG
	frameIterOffset= 0; frameIterMult  = 1;
	frameIterStart = 0; frameIterEnd   = 1;
elif dataSet==201 or dataSet==202 or dataSet==203:
	# 201 = 4 drops low
	# 202 = 6 drops higher
	# 203 = 6 drops higher , 2stepOf
	dat0 = 1; dat1 = 2; # fixed for grids!
	datPath = "pcgDataGrid01/phito02" 
	doLoadGrid = 1;
	dim        = 3;
	resLoad    = 80;
	res        = 96;
	frameIterOffset= 0; frameIterMult  = 1;
	frameIterStart = 0; frameIterEnd   = 1;

	#res        = 34; # small! 
	res        = 44; 
	#res        = 54; 
	#res        = 96; 

	# modify load
	if dataSet==202:
		dataSetSuffix = "B_"
	# spcecial = 2 step OF!
	if dataSet==203:
		#dataSetSuffix = "C_" # org, fastMarch = False
		dataSetSuffix = "D_"  # t2, fastMarch = True


else:
	print( "Error , no valid dataSet chosen!"); exit(1);

# automatically compute border offset
loadOffset  = Vec3(0.)
if autoBorder>0:
	resOffset   = int(res * autoBorder);
	loadOffset  = Vec3(resOffset,resOffset,resOffset);

# repeat the first frame?
doRepeat = False; 
repeatLen = int(res*autoBorder*1.0); 
if 1 and (dataSet>=200 or dataSet==102 or dataSet==103 or dataSet==104): #  NT_DEBUG , check repeat for long
	doRepeat = True;
	loadOffset.z = loadOffset.z + repeatLen


# ====== other params ... ====

# maximal region around surface to include from levelset
# modified below...
maxDist = 5

# solving accuracy, default
cgAccuracy = 1e-04

# cfl condition for advection (old 10)
cfl = 999; # NT_DBEUG off 

# dont correct outer layer in levelset 
corrOuter = False;
# use fast-marching , or simple-reinit?
#doFastMarch = True;
doFastMarch = False; # default for now

# misc params, to be initialized
lsFactor = 0
wSmooth = 0
wEnergy = 0
wPostBlur = 0
maxDist = 0
blurSelection = 0

# for projection step (used with doLsProject)
projMaxDist = 4.0; # 4 & 40 - good defaults...
projMaxIter = 40;

gsStd  = Vec3(res,res,res*thirdDimFactor)
if (dim==2):
	gsStd.z=1

solv           = Solver(name='main', gridSize = gsStd, dim=dim )
solv.timestep  = 1.0
times          = Timings()

# setup relics  ...
flags = solv.create(FlagGrid)
flags.initDomain()
flags.setConst(TypeFluid) # set all to TypeFluid

# testing , preview!
vDiff = solv.create(VecGrid)
i0blend  = solv.create(LevelsetGrid)
i0reinit = solv.create(LevelsetGrid)

# allocate grid data
i0      = solv.create(LevelsetGrid)
i1      = solv.create(LevelsetGrid)
i0org   = solv.create(LevelsetGrid)
i1org   = solv.create(LevelsetGrid)
vel     = solv.create(VecGrid)
velTmp  = solv.create(VecGrid);
# debug multi level solve
debvel5 = solv.create(VecGrid)
debvel4 = solv.create(VecGrid)
debvel3 = solv.create(VecGrid)
debvel2 = solv.create(VecGrid)
debvel1 = solv.create(VecGrid)
debls5 = solv.create(LevelsetGrid)
debls4 = solv.create(LevelsetGrid)
debls3 = solv.create(LevelsetGrid)
debls2 = solv.create(LevelsetGrid)
debls1 = solv.create(LevelsetGrid)
#rhs   = solv.create(RealGrid)
deberr1 = Vec3(0,0,0)
deberr2 = Vec3(0,0,0)
deberr3 = Vec3(0,0,0)
deberr4 = Vec3(0,0,0)
deberr5 = Vec3(0,0,0)
deberr6 = Vec3(0,0,0)
deberr7 = Vec3(0,0,0)
deberr8 = Vec3(0,0,0)
# optional for 2-param interpol
velThird   = 0
iDiffThird = 0

# only for velocity extrapol test , shouldnt be needed in general
velMac = ""

# test motion
i0adv   = solv.create(LevelsetGrid) # advected

# multi step OF
velTmp2  = solv.create(VecGrid)
iTmp1    = solv.create(LevelsetGrid)

iDiff     = solv.create(RealGrid)
iDiffBack = solv.create(RealGrid)

# other solver data
meshDest = solv.create(Mesh)
meshRef1 = solv.create(Mesh)
meshRef0 = solv.create(Mesh)

# solver object for load & size convert
gsload = vec3(resLoad,resLoad,resLoad*thirdDimFactor)
if (dim==2):
	gsload.z=1
sload  = Solver(name='loaddata', gridSize = gsload, dim=dim)
iload  = sload.create(RealGrid)

# 3d mesh env for 2d data
doMeshFrom2d = True
if doMeshFrom2d:
	t2gs = vec3(res,res,10)
	#if (dim==2): t2gs.z=1
	t2dsolv  = Solver(name='mesh3d', gridSize = t2gs, dim=3)
	t2dphi   = t2dsolv.create(LevelsetGrid)

# create preview mesh for 2d data
def previewMesh(phi, meshOut):
	for i in range(0,10):
		placeGrid2d(phi,t2dphi,dstz=i)
	t2dphi.setBound(0.1,0);
	t2dphi.createMesh(meshOut);

#///

if showUi and (GUI):
	gui = Gui()
	gui.windowSize(dispWidth,dispHeight);
	gui.setCamPos (0,0,-camPosZ);
	gui.show()
	if not regen2dSlices:
		gui.pause();
	else:
		gui.toggleHideGrids();

avgErrBest  = 1e10
screenNum   = 0
blendScrNum = 0


# helper functions
def extrapolLs(target, maxt=-1):
	# extrapolation distance, default & override
	d = (maxDist+1);
	if maxt>=0:
		d = maxt;

	if 1: # for testing, disable extrapol
		if doFastMarch:
			target.reinitMarching(flags=flags, maxTime=d , correctOuterLayer=corrOuter ); 
		else:
			extrapolateLsSimple(phi=target, distance=d)  
			extrapolateLsSimple(phi=target, distance=d, inside=true)

	if doConstOverride: # for testing, turn into 0-1 surface
		print("WARNING - const override test active!");
		constOverride(phi=target, val=-(1./lsFactor) );

	# extrapolLs done

# process newly loaded data for of
def doPostProcess(iTarget):
	if(isLevelset):
		extrapolLs(iTarget);
		remapLsValues( grid=iTarget, max=maxDist ) 
		iTarget.multConst(lsFactor)
	else:
		# convert to LS:
		iTarget.multConst(-1.0)
		iTarget.addConst(0.5) # org
		iTarget.addConst(0.3) # larger off

		# workaround , convert to LS...
		if 1:
			extrapolLs(iTarget);
			remapLsValues( grid=iTarget, max=maxDist )
			iTarget.multConst(lsFactor)
		#unused: simpleBlur(iTarget, preImgBlur)

	#iTarget.multConst(-1.0); 
	#iTarget.setBound( int(maxDist*-0.25), int(res*autoBorder*0.5) ); # old, hard reset
	iTarget.setBoundNeumann(int(res*autoBorder*0.5));
	# end post process
	

# load and process data file
def getDataName(dataId, num, attachSuffix=True ):
	suffix = "";

	# hard coded offset for paper 2d example
	# try1
	if 0:
		if dataSet==6 and dataId==1:
			num = 100 # img 7
		if dataSet==6 and dataId==2:
			num = 70 + 10
	# try2
	if 1:
		if dataSet==6 and dataId==1: 
			#num = 60 
			num = num + 20 
		if dataSet==6 and dataId==2:
			#num = 40 
			num = num + 0 
	# end - hard coded offset for paper 2d example

	if attachSuffix:
		suffix = ".uni";
	if(doLoadGrid<=0):
		return("./%s%02d_%04d%s" % (datPath, dataId, num, suffix ) );
	else:
		# grid data
		if  (dataId==1):
			return("./%s_x%dy%dz%d_%04d%s" % (datPath, gridLoadX,gridLoadY,gridLoadZ ,num, suffix ) );
		elif(dataId==2):
			return("./%s_x%dy%dz%d_%04d%s" % (datPath, gridOtherX,gridOtherY,gridOtherZ ,num, suffix ) );
		elif(dataId==3):
			return("./%s_x%dy%dz%d_%04d%s" % (datPath, gridThirdX,gridThirdY,gridThirdZ ,num, suffix ) );
		else:
			print( "Error - invalid grid dataId %d \n"%(dataId) ); 
			exit(1);
	return("-error-");

def loadData(iTarget, dataId, num ):
	global iload;
	if 1:
		# regular load
		iload.load( getDataName(dataId, num) );

		if dataSet==6: # 2d paper example
			iload.setBound(0.1, 3 ); # paper example - cut off flat side

	else:
		# fixed init of sphere for debugging
		drop = solv.create(Sphere, center=vec3(0.4,0.5,0.0)*vec3(res,res,res), radius=res*0.25); 
		if dataId==2:
			drop = solv.create(Sphere, center=vec3(0.4+0.05,0.5,0.0)*vec3(res,res,res), radius=res*0.25); 
		iload = drop.computeLevelset()
		# fixed init of sphere, end

	interpolateGrid(iTarget, iload, offset=loadOffset , scale=vec3(1.-2*autoBorder) );

	if doRepeat: 
		repeatFrame( iTarget, loadOffset.z + 1./float(res), repeatLen, bnd=0 ); 

	iTarget.setBound(0.1, int(res*autoBorder) ); # reset border of LS! 


def getResultName(dataId0, dataId1, num, suffix ):
	name = "";
	if(doLoadGrid<=0):
		name = "%s_res%02d_%s.uni" % (getDataName(dataId0, num, False), dataId1, suffix);
	else:
		name = "%s_to_x%dy%dz%d_%s.uni" % (getDataName(dataId0, num, False), gridOtherX,gridOtherY,gridOtherZ, suffix);
	return name;


def storeResult(iTarget, dataId0, dataId1, num , suffix):
	mySuffix = dataSetSuffix+suffix
	if doStoreResult:
		iTarget.save( getResultName(dataId0, dataId1, num, mySuffix) );

def loadResult(iTarget, dataId0, dataId1, num , suffix):
	mySuffix = dataSetSuffix+suffix
	iTarget.load( getResultName(dataId0, dataId1, num, mySuffix) );


# small helper to prepare final LS
def createFinal(iTarget, iSrc, copyFromSrc=True ):
	if(copyFromSrc):
		iTarget.copyFrom(iSrc);
	iTarget.setBound(-0.1, int(res*autoBorder) ); # reset border of LS!
	if(isLevelset):
		dummy = 1.;
		#extrapolLs(iTarget); # only slight gain in accuracy...

# simple helper to select 1 run from grid
def checkSingleSelect(testRunX,testRunY,testRunZ, otherRunX,otherRunY,otherRunZ):
	if(testRunX == singleSel[0] and
	   testRunY == singleSel[1] and
	   testRunZ == singleSel[2] and
	   otherRunX== singleSel[3] and
	   otherRunY== singleSel[4] and
	   otherRunZ== singleSel[5]): 
		return True;
	return False;


############################################################################################################################

accuAvgErr = 0.; # global results of runseq
def runSeq():
	global screenNum, blendScrNum, avgErrBest
	global accuAvgErr
	global gridOtherX, gridOtherY, gridOtherZ, gridThirdX, gridThirdY, gridThirdZ;

	accuAvgErr = 0.;
	accuAvgCnt = 0.;

	#main loop, check different tests
	for frameIter in range(frameIterStart,frameIterEnd): 

		# =========================================== setup

		i0.setConst(0.);
		i1.setConst(0.);
		print( "\nframeIter: %dd "%(frameIter) )

		# load & backup input data
		fileNum = frameIterOffset + frameIter*frameIterMult
		loadData( i0, dat0, fileNum );
		loadData( i1, dat1, fileNum );
		i0org.copyFrom(i0); 
		i1org.copyFrom(i1); 
		if 1 and isLevelset: # why extrapolate original?  -> data should match
			extrapolLs(i0org);
			extrapolLs(i1org);
		doPostProcess(i0);
		doPostProcess(i1);
		
		if not doOnlyLoad:
			# =========================================== step

			vel.setConst(vec3(0,0,0));

			opticalFlowMultiscale3d(vel=vel, i0=i0, i1=i1,   wSmooth=wSmooth, wEnergy=wEnergy  , postVelBlur=wPostBlur 
				, cgAccuracy=cgAccuracy , blurType=blurSelection , cfl=cfl, orderTime=orderTime,orderSpace=orderSpace  
				, resetBndWidth=autoBorder , multiStep=multiStep , projSizeThresh=projSizeThreshold, minGridSize=6, doFinalProject=doFinalProject 
				, dv1=debvel1,dv2=debvel2,dv3=debvel3,dv4=debvel4,dv5=debvel5 , dr1=debls1,dr2=debls2,dr3=debls3,dr4=debls4,dr5=debls5 );

			# use level set projection
			if doLsProject: 
				maxLSPNorm = 0.50; # zero defo where normals disagree (old: skip vary normals & extrapolate if <1 )
				maxLSPDist = 20;   # skip when defo further than this
				lspBlur    = wPostBlur * 0.5;
				#lspBlur    = 0.01; maxLSPDist = 5; # settings, test blur fill in

				# control by loop
				maxLSPNorm = -1. # doesnt matter
				maxLSPDist = projMaxDist;
				#lspBlur    = wPostBlur * 0.5; 
				lspBlur    = 2.
				maxIter    = projMaxIter;

				print( "\nAdding level-set projection, params %d %f , blur %d, iter %d " % (maxLSPDist, maxLSPNorm, lspBlur , maxIter) );

				corrVelsOf3d( velTmp , vel, i0, i0, i1 , maxLSPDist, maxLSPNorm , lspBlur, autoBorder , maxIter ); 
				#vel.addScaled(velTmp, Vec3(-1.0) );

			#if 0:
				#maxVel = vel.getMaxValue(); print( "Max vel "+str(maxVel) )# check... 


			# =========================================== advect forward 
			i0adv.copyFrom(i0);  # calcErrOrg
			advectSemiLagrangeCfl(cfl=cfl, flags=flags, vel=vel, grid=i0adv, order=orderTime)
			i0adv.setBoundNeumann(0); 

			# calcErrOrg 
			createFinal( i0reinit, i0adv );

			# =========================================== calc error
			# load reference surface again
	
			# just for comparison
			if (dim==3) and showUi:
				i1org.createMesh(meshRef1)
				i0org.createMesh(meshRef0)

			avgErr1 = 0.

			# create full SDF for advected i0
			if(isLevelset): 
				# mark disagreements in levelsets (also compute error
				avgErr1 = calcLsDiff3d(i0=i0adv, i1=i1, correction=20. ) * 1.; # calcErrOrg
				#? i0reinit.multConst(-1.); createFinal(i0reinit, i0reinit, False); # calcErrOrg

				if 1 and (dim==3) and showUi:
					iTmp1.copyFrom(i0reinit); 
					iTmp1.multConst( -( maxDist/(lsFactor) ));
					iTmp1.createMesh(meshDest)

			else:
				# non levelsets...
				avgErr1 = calcDensDiff(i0=i0reinit, i1=iTmp1); 

			print( "Frame %d %d->%d: avg curr error %f " % (fileNum, dat0,dat1, avgErr1 ) ) # per frame error

			# =========================================== 


			# recons 3d mesh surface for 2d
			if doMeshFrom2d and dim==2:
				iTmp1.copyFrom(i0reinit); 
				iTmp1.multConst( -( maxDist/(lsFactor) ));
				#simpleBlur(iTmp1, 1.);
				previewMesh(iTmp1, meshDest);
				previewMesh(i1org, meshRef1);
				previewMesh(i0org, meshRef0);


			# store difference grid
			iDiff.copyFrom(i1org); 
			iDiff.sub(i0reinit);
			#iDiff.setBound(0., int(res*autoBorder) ); # reset border 

			# store result
			storeResult( iDiff, dat0,dat1, fileNum, "dif" );
			storeResult( vel,   dat0,dat1, fileNum, "vel" );

		else:
			# load only
			loadResult( iDiff, dat0,dat1, fileNum, "dif" );
			loadResult( vel  , dat0,dat1, fileNum, "vel" );
			avgErr1 = 0.;


		# load data for third 
		if doLoadThird:
			# a bit ugly - bend "other" selector towards third
			(tmpx, tmpy, tmpz) = (gridOtherX, gridOtherY, gridOtherZ)
			(gridOtherX, gridOtherY, gridOtherZ) = getGridDataId(gridThirdX, gridThirdY, gridThirdZ);

			loadResult( iDiffThird, dat0,dat1, fileNum, "dif" );
			loadResult( velThird  , dat0,dat1, fileNum, "vel" );

			(gridOtherX, gridOtherY, gridOtherZ) = (tmpx, tmpy, tmpz)

		accuAvgErr += avgErr1 
		accuAvgCnt += 1.0

		solv.step()
		if writePngs and showUi: # output
			gui.screenshot( 'out02ms3_%05d.png' % screenNum ); 
			screenNum += 1

		# =========================================== 
		# try blending...
		maxBl = blendIters
		for blendIt in range(0,maxBl): 
			if blendSingleAlpha>=0.:
				alpha = blendSingleAlpha;
			else:
				alpha = (1.0*blendIt) / (1.0*(maxBl-1.0));
			beta = blendThirdAlpha; # optional

			i0blend.copyFrom(i0org); 

			# blending: 0 off, 1 SL, 2 wDiff , 3 2param
			blendMode = 1
			if not doLoadThird:
				print( "Blending %f "%(alpha) )
			else:
				print( "Blending sec %f, third %f "%(alpha, beta) )
				blendMode = 3;

			# show org i1 last
			if alpha>(1.0+1e-04): 
				blendMode = 99;

			# ... blend!
			if blendMode==0:
				# naive LS blend
				i0blend.multConst( 1.0 - alpha )
				iTmp1.copyFrom(i1org);
				iTmp1.multConst( alpha )
				i0blend.add( iTmp1 );
			elif blendMode==1:
				# just forward
				advectSemiLagrangeCfl(cfl=cfl, flags=flags, vel=vel, grid=i0blend, order=orderTime, velFactor=alpha)
			elif blendMode==2:
				# forward & backw
				advectSemiLagrangeCfl(cfl=cfl, flags=flags, vel=vel, grid=i0blend, order=orderTime, velFactor=alpha )
				if 1 and (isLevelset):
					createFinal( i0blend, 0 , copyFrom_i0adv=False ); # damn, needed?
				diffDuration = 0.3
				#wdiff = 1.0 * alpha # org identical
				wdiff = (1.0/diffDuration) * alpha - (1.0/diffDuration-1.) # diffDuration% shift
				# only add if >0
				if wdiff>0.:
					iTmp1.copyFrom( iDiff );   
					iTmp1.multConst( wdiff )
					advectSemiLagrangeCfl(cfl=cfl, flags=flags, vel=vel, grid=iTmp1, order=orderTime, velFactor=-(1.-alpha) )
					i0blend.add( iTmp1 );
					print( "Adding diff vector, weight = %f "%(wdiff) )

			elif blendMode==3:
				# 2-param, just forward
				advectSemiLagrangeCfl(cfl=cfl, flags=flags, vel=vel, grid=i0blend, order=orderTime, velFactor=alpha)
				advectSemiLagrangeCfl(cfl=cfl, flags=flags, vel=velThird, grid=i0blend, order=orderTime, velFactor=beta)

			elif blendMode==99:
				print( "No blend! Showing org1");
				i0blend.copyFrom(i1org); 

			else:
				print( "Invalid blend mode! %d"%blendMode ); exit(1);



			# reinit...
			if(isLevelset):
				extrapolLs(i0blend, res/4);
				i0blend.setBound(0.1, int(res*autoBorder) ); # reset border 
				if (dim==3) and showUi:
					i0blend.createMesh(meshDest)
				if (dim==2) and showUi:
					previewMesh(i0blend, meshDest);
				# optionally save mesh...?
				if 0: # output
					meshDest.save( 'mesh_%d_%04d.obj' % (blendIt,blendScrNum) )
			else:
				i0blend.copyFrom(i0blend, False)
				if 0: # output
					i0blend.save( 'outblxy_%d_%04d.uni' % (blendIt,blendScrNum) )

			# output


			# re-generate 2d slices for whole 2d anim, output - similar to preview Mesh , only for solves with time dim
			if regen2dSlices:
				# note - skip region that is only safety border (ie, res* autoBorder*(1+1) )
				# once for blank region, and once in the beginning for the first-frame-repeat
				for z in range( int(res* autoBorder*2.0)-1, int(res*thirdDimFactor-(res* autoBorder)) ):
					getSliceFrom3d(src=i0blend, srcz=z, dst=t2dphi, dstz=0);
					for dstz in range(1,10):
						getSliceFrom3d(src=t2dphi, srcz=0, dst=t2dphi, dstz=dstz);
					t2dphi.setBound(0.1,0);
					t2dphi.createMesh(meshDest);
					if writePngs and showUi: # output
						#gui.screenshot( 'out_%d_%04d.png' % (blendIt,z) ); 
						if not doLoadThird:
							gui.screenshot( 'outbl%02d_%f_%04d.png' % (dataSet,alpha,z) ); 
						else:
							gui.screenshot( 'outbl%02d_%f_%f_%04d.png' % (dataSet,alpha,beta,z) ); 
					solv.step();
			else:
				# normal display / screenshot
				solv.step()
				if 0 and writePngs and showUi: # output
					gui.screenshot( 'outbl01_%d_%04d.png' % (blendIt,blendScrNum) ); 


		if maxBl>0:
			print("Blending done");
			solv.step(); # stop after blending is done?
		# done blending 

		blendScrNum = blendScrNum + 1


	accuAvgErr /= accuAvgCnt
	print( "Accu avg error %f  (params lsf %f, maxd %f, smo %f, ene %f, blr %f, cgacc %f, rI %d, rO %d)   d/pSet:%d,%d|%d,%d " % (accuAvgErr, lsFactor, maxDist, wSmooth, wEnergy, wPostBlur, cgAccuracy, res, resLoad, dataSet, paramSet, dat0,dat1 ) )
	if(accuAvgErr < avgErrBest):
		avgErrBest = accuAvgErr
	times.display()
	exit(1); # quit after one cycle 
		

############################################################################################################################

def loopParamSearch() :
	global lsFactor, wSmooth, wEnergy, wPostBlur;
	global maxDist, cgAccuracy, blurSelection;
	global projMaxDist , projMaxIter;

	# prepare loop ranges
	loopLenMd  = 1
	loopLenLfs = 1
	loopLenBlr = 1

	if paramSet==6: # search
		loopLenMd  = 3
		loopLenLfs = 5
		loopLenBlr = 5

	if paramSet==7: # search
		loopLenMd  = 5
		loopLenLfs = 5
		loopLenBlr = 5


	# main loop, vary parameters
	for testRunMd in range(0,loopLenMd):
		for testRunLsf in range(0,loopLenLfs): # 5
			for testRunBlr in range(0,loopLenBlr): # 10

				# run1, main params (range 4,6,6)
				if paramSet==1: # search 64
					lsFactor  = pow(10, (1.-testRunLsf) )
					wSmooth   = 0.00002
					wEnergy   = 0.001
					wPostBlur = 1.0 * testRunBlr
					maxDist   = 0.5 * pow(5, (testRunMd) ) 

				# good run 102 accu avg error 3.562538  (params lsf 0.100000, maxd 12.500000, smo 0.000020, ene 0.001000, blr 3.000000)
				elif paramSet==2: # search 64 p2
					lsFactor  = 0.1
					wSmooth   = 0.00002
					wEnergy   = 0.001
					wPostBlur = 3.0 
					maxDist   = 12.0
					# run2
					wSmooth   = 0.000001 * pow(5, testRunBlr)
					wEnergy   = 0.0001   * pow(5, testRunLsf)

				elif paramSet==3: # search 140
					lsFactor  = 0.1
					wEnergy   = 0.000125
					wPostBlur = 3.0 * testRunBlr
					maxDist   = 0.5 * pow(5, (testRunMd) ) 
					wSmooth   = 0.000001 * pow(5, testRunLsf)

				# 1/25 of smoothing amount, 3x more blur!
				elif paramSet==4:
					lsFactor  = 0.1
					wSmooth   = 0.000001
					wEnergy   = 0.000125
					wPostBlur = 9.0 
					maxDist   = 12.0

				# good from 64 , smoke
				elif paramSet==5:
					wPostBlur  = 3.0 
					wSmooth    = 0.01
					wEnergy    = 0.0005
					# 1 for median, 2 for simple
					#wPostBlur  = 2.0  

					wPostBlur  = 1.0 * testRunBlr
					blurSelection  = testRunLsf

					#wPostBlur  = 1.0 
					blurSelection  = 1

					maxDist    = 1. + 10. *testRunLsf
					lsFactor   = 0.1 * (testRunBlr+1)
					maxDist    = 20; 
					lsFactor   = 0.1
					wPostBlur  = 6. # compromise?

					#cgAccuracy = 0.003
					cgAccuracy = 0.03 

				# exhaustive search
				elif paramSet==6:
					maxDist    = 20; 
					lsFactor   = 0.1;
					cgAccuracy = 0.01;
					blurSelection = 1;
					wPostBlur = 2.0       *         testRunBlr;
					wSmooth   = 0.00001   * pow(10, testRunLsf);
					wEnergy   = 0.0001    * pow(10, testRunMd);

				# search project
				elif paramSet==7:
					wSmooth    = 0.01
					wEnergy    = 0.0005 
					blurSelection  = 1 
					maxDist    = 20; 
					lsFactor   = 0.1
					wPostBlur  = 6. 
					cgAccuracy = 0.001 
					wPostBlur   = 10. / (testRunBlr+1);
					projMaxDist = 20  / (testRunLsf+1);
					projMaxIter = 200 / (testRunMd+1);
					# search 1: maxdist 20 bad!
					wPostBlur   = 2.* (testRunBlr+0.);
					projMaxDist = 5-(testRunLsf);
					projMaxIter = 5 * (testRunMd+1);
					# good:
					wPostBlur   = 4;  # times 0.5 later -> 2
					projMaxIter = 40;
					projMaxDist = 4;
					# test fill in , NT_DEBUG , doFillIn
					# minimal difference?
					projMaxIter = 3 * (testRunBlr+0);

				# liquid 64 3d
				elif paramset==10:
					wSmooth    = 0.000025
					wEnergy    = 0.000125
					wPostBlur  = 3.0 
					maxDist    = 12.0
					lsFactor   = 0.1
					#cgAccuracy = pow(10, -testRunBlr)
					#cgAccuracy = 0.05
					cgAccuracy = 0.5   # low, but still ok...
					#cgAccuracy = 10.0 # solve effectively off

				else:
					print( "Error , no valid paramset chosen!"); exit(1);

				print( "\nparamSearch loop %d %d %d "%(testRunBlr,testRunMd,testRunLsf) )
				runSeq();

			print( "Accu avg error === " ) # loop sep

	print( "Accu avg error best so far %f " % avgErrBest )


############################################################################################################################

def idString(testRunX,testRunY,testRunZ, otherRunX,otherRunY,otherRunZ):
	ret = "";
	# uni directional graph
	#if   otherRunZ < testRunZ:
		#ret = ("%d_%d_%d vs %d_%d_%d" % (otherRunX,otherRunY,otherRunZ, testRunX,testRunY,testRunZ   ));
	#elif otherRunY < testRunY:
		#ret = ("%d_%d_%d vs %d_%d_%d" % (otherRunX,otherRunY,otherRunZ, testRunX,testRunY,testRunZ   ));
	#elif otherRunX < testRunX:
		#ret = ("%d_%d_%d vs %d_%d_%d" % (otherRunX,otherRunY,otherRunZ, testRunX,testRunY,testRunZ   ));
	#else:
		#ret = ("%d_%d_%d vs %d_%d_%d" % (testRunX,testRunY,testRunZ,    otherRunX,otherRunY,otherRunZ));

	# bi-directional
	ret = ("%d_%d_%d vs %d_%d_%d" % (testRunX,testRunY,testRunZ,    otherRunX,otherRunY,otherRunZ));

	return ret;

def getGridDataId(ix,iy,iz):
	outx,outy,outz = (-1,-1,-1)
	if dataSet==200 or dataSet==201: 
		outx = ix * 2;
		outy = iy * 2;
		outz = iz * 2;
	elif dataSet==202:
		outx = 1 + ix * 2;
		outy = 1 + iy * 2;
		outz = 0; # 2 params
	elif dataSet==203:
		outx = 1 + ix * 2;
		outy = 1 + iy * 2;
		outz = 0; # 2 params
	return (outx,outy,outz)

def loopGrid() :
	global lsFactor, wSmooth, wEnergy, wPostBlur;
	global maxDist, cgAccuracy, blurSelection;
	global gridLoadX, gridLoadY, gridLoadZ;
	global gridOtherX, gridOtherY, gridOtherZ;
	global gridThirdX, gridThirdY, gridThirdZ;
	global doLoadThird, velThird, iDiffThird;

	# sanity check for grid single load
	#if doSingleSelect and not doOnlyLoad:
		#print ("Error, use single load only with load mode"); exit(1);

	errAllRuns = 0.
	# prepare loop ranges
	loopLenZ = 1;
	loopLenY = 1;
	loopLenX = 1;

	if dataSet==200 or dataSet==201: 
		loopLenZ = 1
		loopLenY = 2
		loopLenX = 2
	if dataSet==202:
		loopLenZ = 1
		loopLenY = 3
		loopLenX = 2
	if dataSet==203:
		loopLenZ = 1
		loopLenY = 3
		loopLenX = 2

	gridThirdX, gridThirdY, gridThirdZ = (singleSel[6],singleSel[7],singleSel[8]) 
	if doSingleSelect and gridThirdX>=0 and gridThirdY>=0 and gridThirdZ>=0:
		print("Loading third data point for interpolation")
		doLoadThird = True;
		# allocate grids
		velThird   = solv.create(VecGrid)
		iDiffThird = solv.create(LevelsetGrid)

	# fix params
	maxDist   = 20
	lsFactor  = 0.1
	wSmooth   = 0.01
	wEnergy   = 0.0005
	#wPostBlur = 2
	wPostBlur = 6
	#cgAccuracy = 0.003 # doesnt pay off
	cgAccuracy = 0.03
	blurSelection = 1
	#cgAccuracy = 11.11 # NT DEBUG off

	gv = {};

	# main loop, vary parameters
	for testRunZ in range(0,loopLenZ):
		for testRunY in range(0,loopLenY): 
			for testRunX in range(0,loopLenX): 
				# calculate offsets into data
				gridLoadX, gridLoadY, gridLoadZ = getGridDataId(testRunX,testRunY,testRunZ );

				for otherRunZ in range(0,loopLenZ):
					for otherRunY in range(0,loopLenY): 
						for otherRunX in range(0,loopLenX): 
							if 1 and testRunX==otherRunX and testRunY==otherRunY and testRunZ==otherRunZ:
								continue;
							#gridOtherX = otherRunX * 2; gridOtherY = otherRunY * 2; gridOtherZ = otherRunZ * 2;
							gridOtherX, gridOtherY, gridOtherZ = getGridDataId(otherRunX,otherRunY,otherRunZ );

							#print( "\ngridSearch (%d %d %d) -> (%d %d %d)  "%(testRunX,testRunY,testRunZ, otherRunX,otherRunY,otherRunZ) )
							#print( "    sel    (%d %d %d) -> (%d %d %d)  "%(singleSel[0],singleSel[1],singleSel[2],singleSel[3],singleSel[4],singleSel[5]) );
							if (not doSingleSelect) or (doSingleSelect and checkSingleSelect(testRunX,testRunY,testRunZ, otherRunX,otherRunY,otherRunZ) ):
								#print( "\ngridSearch (%d %d %d) -> (%d %d %d)  "%(testRunX,testRunY,testRunZ, otherRunX,otherRunY,otherRunZ) )
								print( "\ngridSearch (%d %d %d) -> (%d %d %d)  "%(gridLoadX,gridLoadY,gridLoadZ, gridOtherX,gridOtherY,gridOtherZ) )

								runSeq();

								# store error
								print( idString(testRunX,testRunY,testRunZ, otherRunX,otherRunY,otherRunZ) );
								gv[ idString(testRunX,testRunY,testRunZ, otherRunX,otherRunY,otherRunZ) ] = accuAvgErr;
								errAllRuns = errAllRuns + accuAvgErr;
				# sub loop done
	# main loop done

	if not doOnlyLoad:
		if 0:
			print( "\nResults:");
			for item in gv:
				print("%s: %f" % (item, gv[item]));

		if 1:
			sorted_list = [x for x in gv.iteritems()] 
			sorted_list.sort(key=lambda x: x[0]) # sort by id string
			#sorted_list.sort(key=lambda x: x[1]) # sort by value 
			#sorted_list.reverse() # to reverse the sort
			print( "\nResults sorted:");
			for item in sorted_list:
				print("%s: %f" % (item[0], item[1]));

	print( "\nTotal added err: %f"%(errAllRuns));
	# grid done



############################################################################################################################


# select loop type
if dataSet<200:
	loopParamSearch();
else:
	loopGrid();