text:
# using  Test, Revise,CUDA
# includet("C:\\GitHub\\GitHub\\NuclearMedEval\\src\\structs\\BasicStructs.jl")

# includet("C:\\GitHub\\GitHub\\NuclearMedEval\\src\\utils\\CUDAGpuUtils.jl")
# includet("C:\\GitHub\\GitHub\\NuclearMedEval\\src\\distanceMetrics\\Housdorff\\mainHouseDorffKernel\\HFUtils.jl")
# includet("C:\\GitHub\\GitHub\\NuclearMedEval\\src\\distanceMetrics\\Housdorff\\mainHouseDorffKernel\\ProcessPadding.jl")
# includet("C:\\GitHub\\GitHub\\NuclearMedEval\\src\\distanceMetrics\\Housdorff\\mainHouseDorffKernel\\ProcessMainData.jl")
# includet("C:\\GitHub\\GitHub\\NuclearMedEval\\test\\GPUtestUtils.jl")

# using ..CUDAGpuUtils, ..HFUtils,Cthulhu

# using ..CUDAGpuUtils, ..HFUtils,Cthulhu

# #@testset "processPaddingTest" begin

#     testArrInCPU= falses(100,100,100);
#     #testArrInCPU[CartesianIndex(1,1,1)]= true;
#     testArrInCPU[CartesianIndex(1,1,1)]= true
#     testArrInCPU[CartesianIndex(5,5,5)]= true
#     testArrIn = CuArray(testArrInCPU);

#     metaDataCPU = falses(15,15,15,4); # all neighbouring set to be inactive and not full
#     metaData= CuArray(metaDataCPU);
#     debugArr= CUDA.zeros(Bool,9);
#     resArraysCounter=CUDA.zeros(Int32,1);

#     referenceArray= CUDA.ones(Bool,100,100,100);
#     resArray = CUDA.zeros(UInt16,50,50,50);
#     blockBeginingX,blockBeginingY,blockBeginingZ =UInt8(33),UInt8(33),UInt8(33);
#     currBlockX,currBlockY,currBlockZ =UInt8(2),UInt8(2),UInt8(2);
#     isPassGold = true;
#     #metaDataCPU[2,2,3,:].= true

#     metadataDims = size(metaData);
#     mainQuesCounter=CUDA.zeros(Int32,1);
#     mainWorkQueue = CUDA.zeros(UInt8,10,4);#simulating place for 10 blocks
#     iterationNumber= UInt32(2);
# function testKernelForPaddingAnalysis(analyzedArr, refAray,blockBeginingX,blockBeginingY,blockBeginingZ,resArray,resArraysCounter
#                                     ,isPassGold,currBlockX,currBlockY,currBlockZ,metaData,metadataDims,mainQuesCounter,mainWorkQueue,iterationNumber,debugArr)
#         resShmem =  @cuStaticSharedMem(Bool,(34,34,34))#+2 in order to get the one padding
#         ProcessMainData.executeDataIterFirstPassWithPadding(analyzedArr, refAray,blockBeginingX,blockBeginingY,blockBeginingZ,resShmem,resArray,resArraysCounter,currBlockX,currBlockY,currBlockZ,isPassGold,metaData,metadataDims,mainQuesCounter,mainWorkQueue,iterationNumber,debugArr)
#         #clearMainShmem(shmem)
#         #now we need to deal with padding in shmem res
#         # processAllPaddingPlanes(blockBeginingX,blockBeginingY,blockBeginingZ,resShmem
#         #                     ,currBlockX,currBlockY,currBlockZ
#         #                     ,analyzedArr,refAray,resArray
#         #                     ,metaData,metadataDims
#         #                     ,isPassGold)

#         return
# end

# @cuda threads=(16,16) blocks=1 testKernelForPaddingAnalysis(testArrIn,referenceArray,blockBeginingX,blockBeginingY,blockBeginingZ,  resArray,
#   resArraysCounter,isPassGold,currBlockX,currBlockY,currBlockZ,metaData,metadataDims,mainQuesCounter,mainWorkQueue,iterationNumber,debugArr)
#   resArraysCounter[1]
#   Int64(maximum(mainWorkQueue))

#   indicies = CartesianIndices(mainWorkQueue);
#   filtered = filter(ind-> mainWorkQueue[ind]>0,indicies  )

# debugArr[5]
# Int64(mainWorkQueue[2,4])

# Int64(maximum(mainWorkQueue))

# indicies = CartesianIndices(mainWorkQueue);
# filtered = filter(ind-> mainWorkQueue[ind]>0,indicies  )

# metaData[2, 2, 3, 2]==true

# debugArr[1]
# debugArr[2]
# debugArr[3]
# debugArr[4]
# debugArr[5]
# debugArr[6]

# testArrInCPU= falses(100,100,100);
# #testArrInCPU[CartesianIndex(1,1,1)]= true;
# testArrInCPU[CartesianIndex(5,5,5)]= true;
# testArrInCPU[5,5,32]=true

# testArrIn = CuArray(testArrInCPU);
# referenceArray= CUDA.ones(Bool,100,100,100);
# resArray = CUDA.zeros(UInt32,32,32,32);
# blockBeginingX,blockBeginingY,blockBeginingZ =UInt8(0),UInt8(0),UInt8(0);
# currBlockX,currBlockY,currBlockZ =UInt8(2),UInt8(2),UInt8(2);
# isPassGold = true;
# metaDataCPU = falses(15,15,15,4); # all neighbouring set to be inactive and not full
# #metaDataCPU[2,2,3,:].= true

# metaData= CuArray(metaDataCPU);
# metadataDims = size(metaData);
# mainQuesCounter=CUDA.zeros(Int32,1);
# mainWorkQueue = CUDA.zeros(Bool,10,4);#simulating place for 10 blocks
# iterationNumber= UInt32(2);

# @device_code_warntype interactive=true @cuda testKernelForPaddingAnalysis(testArrIn,referenceArray,blockBeginingX,blockBeginingY,blockBeginingZ,  resArray,resArraysCounter,isPassGold,currBlockX,currBlockY,currBlockZ,metaData,metadataDims,mainQuesCounter,mainWorkQueue,iterationNumber,debugArr)

# #end

###########loadMainValues
using Revise, Parameters, Logging, Test
using CUDA
includet("C:\\GitHub\\GitHub\\NuclearMedEval\\test\\includeAllUseFullForTest.jl")
using ..CUDAGpuUtils ,..IterationUtils,..ReductionUtils , ..MemoryUtils,..CUDAAtomicUtils
using ..BitWiseUtils,..ResultListUtils, ..MetadataAnalyzePass,..MetaDataUtils,..WorkQueueUtils,..ProcessMainDataVerB,..HFUtils, ..ScanForDuplicates

mainArrGPU= CUDA.zeros(UInt32, 50,50,20)

dataBdim= (32,10,32)

xMeta,yMeta,zMeta = 1,1,1

resShmemblockData= CUDA.zeros(UInt32, 32,10);
shmemblockData= CUDA.zeros(UInt32, 32,10);
shmemPaddings= CUDA.zeros(Bool, 32,32,6);

rowOne = 0
@setBitTo(rowOne,1,true)
@setBitTo(rowOne,5,true)
@setBitTo(rowOne,32,true)

mainArrGPU[1,1,1]= rowOne
mainArrGPU[1,10,1]= rowOne
mainArrGPU[32,1,1]= rowOne

function testLoadMainValues(shmemPaddings,shmemblockData,resShmemblockData,dataBdim,mainArrGPU,xMeta,yMeta,zMeta)

  sync_threads()

  @loadMainValues(mainArrGPU,xMeta,yMeta,zMeta)

    return
end

@cuda threads=(32,10) blocks=1 testLoadMainValues(shmemPaddings,shmemblockData,resShmemblockData,dataBdim,mainArrGPU,xMeta,yMeta,zMeta)

@test shmemblockData[1,1]==rowOne
@test shmemblockData[1,10]==rowOne
@test shmemblockData[32,1]==rowOne

@test isBit1AtPos(resShmemblockData[1,1],1)
@test isBit1AtPos(resShmemblockData[1,1],2)
@test !isBit1AtPos(resShmemblockData[1,1],3)
@test isBit1AtPos(resShmemblockData[1,1],4)
@test isBit1AtPos(resShmemblockData[1,1],5)
@test isBit1AtPos(resShmemblockData[1,1],6)
@test isBit1AtPos(resShmemblockData[1,1],31)
@test isBit1AtPos(resShmemblockData[1,1],32)

@test isBit1AtPos(resShmemblockData[1,10],1)
@test isBit1AtPos(resShmemblockData[1,10],2)
@test isBit1AtPos(resShmemblockData[1,10],4)
@test isBit1AtPos(resShmemblockData[1,10],5)
@test isBit1AtPos(resShmemblockData[1,10],6)
@test isBit1AtPos(resShmemblockData[1,10],31)
@test isBit1AtPos(resShmemblockData[1,10],32)

@test isBit1AtPos(resShmemblockData[32,1],1)
@test isBit1AtPos(resShmemblockData[32,1],2)
@test isBit1AtPos(resShmemblockData[32,1],4)
@test isBit1AtPos(resShmemblockData[32,1],5)
@test isBit1AtPos(resShmemblockData[32,1],6)
@test isBit1AtPos(resShmemblockData[32,1],31)
@test isBit1AtPos(resShmemblockData[32,1],32)

# holding data about result top, bottom, left right , anterior, posterior ,  paddings

@test shmemPaddings[1,1,1]
@test shmemPaddings[1,10,1]
@test shmemPaddings[32,1,1]

@test shmemPaddings[1,1,2]
@test shmemPaddings[1,10,2]
@test shmemPaddings[32,1,2]

@test shmemPaddings[1,1,3]
@test shmemPaddings[5,1,3]
@test shmemPaddings[32,1,3]

@test shmemPaddings[1,10,3]
@test shmemPaddings[5,10,3]
@test shmemPaddings[32,10,3]

@test shmemPaddings[1,1,4]
@test shmemPaddings[5,1,4]
@test shmemPaddings[32,1,4]
@test !shmemPaddings[32,2,4]

@test shmemPaddings[1,1,5]
@test shmemPaddings[1,5,5]
@test shmemPaddings[1,32,5]

@test shmemPaddings[1,1,6]
@test shmemPaddings[1,5,6]
@test shmemPaddings[1,32,6]

#############  validateData
# using Revise, Parameters, Logging, Test
# using CUDA
# includet("C:\\GitHub\\GitHub\\NuclearMedEval\\test\\includeAllUseFullForTest.jl")
# using ..CUDAGpuUtils ,..IterationUtils,..ReductionUtils , ..MemoryUtils,..CUDAAtomicUtils
# using ..BitWiseUtils,..ResultListUtils, ..MetadataAnalyzePass,..MetaDataUtils,..WorkQueueUtils,..ProcessMainDataVerB,..HFUtils, ..ScanForDuplicates

# mainArr= CUDA.zeros(UInt32, 50,50,10)
# refArr= CUDA.zeros(UInt32, 50,50,10)
# targetArr= CUDA.zeros(UInt32, 50,50,10)

# dataBdim= (32,10,32)

# xMeta,yMeta,zMeta = 1,1,1

# resShmemblockData= CUDA.zeros(UInt32, 32,10);
# shmemblockData= CUDA.zeros(UInt32, 32,10);
# shmemPaddings= CUDA.zeros(Bool, 32,32,6);

# threads=(32,10)
# rowOne = 0
# @setBitTo(rowOne,1,true)
# @setBitTo(rowOne,5,true)
# @setBitTo(rowOne,32,true)

# mainArr[1,1,1]= rowOne
# mainArr[1,10,1]= rowOne
# mainArr[32,1,1]= rowOne

# rowB = 0
# @setBitTo(rowB,1,true)
# @setBitTo(rowB,2,true)
# @setBitTo(rowB,5,true)
# @setBitTo(rowB,6,true)
# @setBitTo(rowB,32,true)

# targetArr[1,10,1]= rowB

# blocks =1
# mainArrDims= (50,50,320)
# metaData = MetaDataUtils.allocateMetadata(mainArrDims,dataBdim)
# #metaData = view(MetaDataUtils.allocateMetadata(mainArrDims,dataBdim),1:9,2:3,4:6,: );
# metaDataDims=size(metaData)
# metaSize= metaDataDims
# workQueue,workQueueCounter= WorkQueueUtils.allocateWorkQueue( max(length(metaData),1) )
# metaData[2,2,2,2]=UInt32(1)
# #setting offsets in metadata
# for i in 1:14
#   metaData[xMeta,yMeta,zMeta,getResOffsetsBeg()+i]=i*10
# end
# isGold = 1
# iterNumb = 1
# resList = allocateResultLists(1000,1000)
# paddingStore = CUDA.zeros(UInt8, metaSize[1],metaSize[2],metaSize[3],32,32)
# inBlockLoopXZIterWithPadding = cld(32,10)
# numberToLooFor = 2
# function testvalidateData(numberToLooFor,inBlockLoopXZIterWithPadding,paddingStore,resList,shmemPaddings,shmemblockData,resShmemblockData,metaData,metaDataDims,mainArrDims,isGold,xMeta,yMeta,zMeta,iterNumb,mainArr,refArr,targetArr,dataBdim,workQueue,workQueueCounter)
#   isMaskFull = true
#   @loadMainValues(mainArr,xMeta,yMeta,zMeta)

#   sync_threads()
#   @validateData(isGold,xMeta,yMeta,zMeta,iterNumb,mainArr,refArr,targetArr)

#     return
# end

# @cuda threads=threads blocks=blocks testvalidateData(inBlockLoopXZIterWithPadding,paddingStore,resList,shmemPaddings,shmemblockData,resShmemblockData,metaData,metaDataDims,mainArrDims,isGold,xMeta,yMeta,zMeta,iterNumb,mainArr,refArr,targetArr,dataBdim,workQueue,workQueueCounter)

# @test shmemPaddings[1,1,4]
# @test shmemPaddings[5,1,4]
# @test shmemPaddings[32,1,4]
# @test !shmemPaddings[32,2,4]

# @test !shmemPaddings[32,2,4]

# @setBitTo(rowB,1,true)
# @setBitTo(rowB,2,true)
# @setBitTo(rowB,5,true)
# @setBitTo(rowB,6,true)
# @setBitTo(rowB,32,true)

# targetArr[1,10,1]= rowB

################## execute Data Iter With Padding  first data

using Revise, Parameters, Logging, Test
using CUDA
includet("C:\\GitHub\\GitHub\\NuclearMedEval\\test\\includeAllUseFullForTest.jl")
using ..CUDAGpuUtils ,..IterationUtils,..ReductionUtils , ..MemoryUtils,..CUDAAtomicUtils
using ..Housdorff,..ProcessMainDataVerB,..MainLoopKernel,..BitWiseUtils,..ResultListUtils, ..MetadataAnalyzePass,..MetaDataUtils,..WorkQueueUtils,..ProcessMainDataVerB,..HFUtils, ..ScanForDuplicates

mainArr= CUDA.zeros(UInt32, 300,300,10)
refArr= CUDA.zeros(UInt8, 300,300,320)
numbToLooFor = 2
dataBdim= (32,5,32)

############ some arbitrary data  meta 2,2,2 with some results to be set
xMeta,yMeta,zMeta = 2,2,2

threads=(dataBdim[1],dataBdim[2])
rowOne = UInt32(0)

#below we will organize it so we should have results in all neghbouring blocks - also from corners
@setBitTo(rowOne,1,true)
@setBitTo(rowOne,5,true)
@setBitTo(rowOne,32,true)

mainArr[33,dataBdim[2]+1,2]= rowOne
mainArr[33,dataBdim[2]*2,2]= rowOne
mainArr[64,dataBdim[2]+1,2]= rowOne
mainArr[64,dataBdim[2]*2,2]= rowOne

mainArr[64,dataBdim[2]+1,2]= rowOne

################   here we will get a block that will become full after the dilatation meta 3,3,3
fullOnesTobecome = UInt32(0)
for bitPos in 1:32
 if(isodd(bitPos))
    @setBitTo(fullOnesTobecome,bitPos,true)
 end
end
for xx in ((2*32)+1):((3*32)), yy in ((2*dataBdim[2])+1):((3*dataBdim[2]))
  mainArr[xx,yy,3]= fullOnesTobecome
end

mainArr[100,16,3]

# for xx in ((3*32)+1):((4*32)), yy in ((3*dataBdim[2])+1):((4*dataBdim[2])), zz in ((3*32)+1):((4*32))
#   refArr[xx,yy,zz]= 2
# end

refArr[(dataBdim[1])+1 ,(dataBdim[2])+1,dataBdim[3]+2]=2 #from top
refArr[(dataBdim[1])+1 ,(dataBdim[2])+1,dataBdim[3]+4]=2 #from bottom

refArr[(dataBdim[1]*2)-1 ,(dataBdim[2])+1, dataBdim[3]+5]=2 #from right
refArr[(dataBdim[1])+2 ,(dataBdim[2])+1,dataBdim[3]+5]=2 #from left

refArr[(dataBdim[1])*2 ,(dataBdim[2]*2)-1,dataBdim[3]+5]=2 #from anterior
refArr[(dataBdim[1])+1 ,(dataBdim[2])+2,dataBdim[3]+5]=2 #from posterior

##single voxels dilatation set ref array
refArr[(dataBdim[1]*4)+4 ,(dataBdim[2]*4)+4,(2-1)*dataBdim[3]+1]=2 #from top
refArr[(dataBdim[1]*4)+2 ,(dataBdim[2]*4)+2,(1-1)*dataBdim[3]+32]=2 #from bottom
refArr[(dataBdim[1]*4) ,(dataBdim[2]*4)+2,(3-1)*dataBdim[3]+5]=2 #from right
refArr[(dataBdim[1]*5)+1 ,(dataBdim[2]*4)+2,(4-1)*dataBdim[3]+5]=2 #from left
refArr[(dataBdim[1]*4)+2 ,(dataBdim[2]*4),(5-1)*dataBdim[3]+5]=2 #from anterior
refArr[(dataBdim[1]*4)+2 ,(dataBdim[2]*5)+1,(6-1)*dataBdim[3]+5]=2 #from posterior

(dataBdim[1]*4)+4
(dataBdim[2]*4)+4
(2-1)*dataBdim[3]+1

(dataBdim[1]*4)+2
(dataBdim[2]*4)+2
1*dataBdim[3]+32
#### single points to check dilatations

#bottom 5,5,1
roww = UInt32(0)
@setBitTo(roww,32,true)
mainArr[(dataBdim[1]*4)+4 ,(dataBdim[2]*4)+4,1]= roww
#top 5,5,2
roww = UInt32(0)
@setBitTo(roww,1,true)
mainArr[(dataBdim[1]*4)+2 ,(dataBdim[2]*4)+2,2]= roww
#left 5,5,3
roww = UInt32(0)
@setBitTo(roww,5,true)
mainArr[(dataBdim[1]*4)+1 ,(dataBdim[2]*4)+2,3]= roww
#right 5,5,4
roww = UInt32(0)
@setBitTo(roww,5,true)
mainArr[(dataBdim[1]*5) ,(dataBdim[2]*4)+2,4]= roww
#posterior 5,5,5
roww = UInt32(0)
@setBitTo(roww,5,true)
mainArr[(dataBdim[1]*4)+2 ,(dataBdim[2]*4)+1,5]= roww
#anterior 5,5,6
roww = UInt32(0)
@setBitTo(roww,5,true)
mainArr[(dataBdim[1]*4)+2 ,(dataBdim[2]*5),6]= roww

### configurations

mainArrDims= size(mainArr)
metaData = MetaDataUtils.allocateMetadata(mainArrDims,dataBdim)
#metaData = view(MetaDataUtils.allocateMetadata(mainArrDims,dataBdim),1:9,2:3,4:6,: );
metaDataDims=size(metaData)

workQueue,workQueueCounter= WorkQueueUtils.allocateWorkQueue( max(length(metaData),1) )
metaData[2,2,2,2]=UInt32(1)
#setting offsets in metadata

twoTwoTwoOffset=0

metaBlock = 0
for xMetaa in 1:6,yMetaa in 1:6, zMetaa in 1:6
  metaBlock+=20
  if(xMetaa==2 && yMetaa==2 && zMetaa==2 )
      twoTwoTwoOffset=metaBlock

  end

  for i in 1:14
    metaData[xMetaa,yMetaa,zMetaa,getResOffsetsBeg()+i]=i*10+metaBlock*15*10
    metaData[(xMetaa),(yMetaa),(zMetaa),(getIsToBeAnalyzedNumb() +i)] =1
  end
end

for xMetaa in 1:6,yMetaa in 1:6, zMetaa in 1:6
  for i in 1:14
    metaData[(xMetaa),(yMetaa),(zMetaa),(getIsToBeAnalyzedNumb() +i)] =1
  end
end

isGold = 1
iterNumb = 1

inBlockLoopXZIterWithPadding = cld(32,dataBdim[2])
numberToLooFor = 2
resList = allocateResultLists(1000000,1000000)
workQueue[:,1] = [2,2,2,1]
workQueue[:,2] = [3,3,3,1]

workQueue[:,3] = [4,3,3,1]
workQueue[:,4] = [2,3,3,1]
workQueue[:,5] = [3,4,3,1]
workQueue[:,6] = [3,2,3,1]
workQueue[:,7] = [3,3,4,1]
workQueue[:,8] = [3,3,2,1]

workQueue[:,9] = [1,2,2,1]
workQueue[:,10] = [3,2,2,1]
workQueue[:,11] = [2,1,2,1]
workQueue[:,12] = [2,3,2,1]
workQueue[:,13] = [2,2,1,1]
workQueue[:,14] = [2,2,3,1]

workQueue[:,15] = [5,5,1,1]
workQueue[:,16] = [5,5,2,1]
workQueue[:,17] = [5,5,3,1]
workQueue[:,18] = [5,5,4,1]
workQueue[:,19] = [5,5,5,1]
workQueue[:,20] = [5,5,6,1]

workQueue[:,21] = [4,5,3,1]
workQueue[:,22] = [6,5,4,1]
workQueue[:,23] = [5,4,5,1]
workQueue[:,24] = [5,6,6,1]

workQueueCounter[1] = 24
dilatationArrs= (mainArr,mainArr)
referenceArrs=(refArr,refArr)
shmemSumLengthMaxDiv4 = 20
paddingStore = CUDA.zeros(UInt8, metaDataDims[1],metaDataDims[2],metaDataDims[3],32,32)
function testProcessDataBlock(numberToLooFor,refArr,inBlockLoopXZIterWithPadding,paddingStore,shmemSumLengthMaxDiv4,referenceArrs,dilatationArrs,resList, metaData,metaDataDims,mainArrDims,isGold,iterNumb,mainArr,dataBdim,workQueue,workQueueCounter)
  MainLoopKernel.@mainLoopKernelAllocations(dataBdim)
  @ifXY 1 1 for i in 1:14
    areToBeValidated[i]=true
  end
  @ifXY 3 1 goldToBeDilatated[1]=true
  @ifXY 3 1 segmToBeDilatated[1]=true

  sync_threads()

  #we get dilatation from block its padding will be analyzed later
  @iterateOverWorkQueue(begin
  ProcessMainDataVerB.@executeDataIter(mainArrDims
      ,dilatationArrs[shmemSum[shmemIndex*4+4]+1]
      ,referenceArrs[shmemSum[shmemIndex*4+4]+1]
      ,shmemSum[shmemIndex*4+1]#xMeta
      ,shmemSum[shmemIndex*4+2]#yMeta
      ,shmemSum[shmemIndex*4+3]#zMeta
      ,shmemSum[shmemIndex*4+4]#isGold
      ,iterationNumberShmem[1]#iterNumb
  )

  end )

  sync_grid(grid_handle)

  @iterateOverWorkQueue(begin
  ProcessMainDataVerB.@executeIterPadding(dilatationArrs[shmemSum[shmemIndex*4+4]+1]
      ,referenceArrs[shmemSum[shmemIndex*4+4]+1]
      ,(shmemSum[shmemIndex*4+1])#xMeta
      ,(shmemSum[shmemIndex*4+2])#yMeta
      ,(shmemSum[shmemIndex*4+3])#zMeta
      ,shmemSum[shmemIndex*4+4]#isGold
      ,iterationNumberShmem[1]#iterNumb
      )

  end )

    return
end

@cuda threads=threads blocks=1 cooperative = true shmem = get_shmemMainKernel(dataBdim) testProcessDataBlock(numberToLooFor,refArr,inBlockLoopXZIterWithPadding,paddingStore,shmemSumLengthMaxDiv4,referenceArrs,dilatationArrs,resList, metaData,metaDataDims,mainArrDims,isGold,iterNumb,mainArr,dataBdim,workQueue,workQueueCounter)

#we need to test couple thing
#1) does dilateted data correctly was written to correct spot in the mainArr
mainArr[1]
referenceArrs[1][1,1,1]
referenceArrs[2][1,1,1]

afterDil = UInt32(0)
for bitPos in [1,2,4,5,6,31,32]
    @setBitTo(afterDil,bitPos,true)
end
afterDil
nn = mainArr[33,dataBdim[2]+1,2]
isBit1AtPos(nn,32)
mainArr[32,dataBdim[2]+1,2]

#bottom 5,5,1
nn = paddingStore[5,5,1,4,4]
@test isBit1AtPos(nn,2)
#top 5,5,2
nn = paddingStore[5,5,2  ,2,2]
@test isBit1AtPos(nn,7)
#left 5,5,3
nn = paddingStore[5,5,3  ,5,2]
@test isBit1AtPos(nn,3)
#right 5,5,4
nn = paddingStore[5,5,4  ,5,2]
@test isBit1AtPos(nn,4)
#posterior 5,5,5
nn = paddingStore[5,5,5  ,2,5]
@test isBit1AtPos(nn,6)
#anterior 5,5,6
nn = paddingStore[5,5,6  ,2,5]
@test isBit1AtPos(nn,5)

###### now those single padding points is it couse dilatations
#bottom - testing top of block below
roww = UInt32(0)
@setBitTo(roww,1,true)
# filter(it->it>0, Int64.(mainArr[(dataBdim[1]*4)+4 ,(dataBdim[2]*4)+4,:]))

@test  mainArr[(dataBdim[1]*4)+4 ,(dataBdim[2]*4)+4,2]== roww
#top - testing bottom of above
roww = UInt32(0)
@setBitTo(roww,32,true)
@test  mainArr[(dataBdim[1]*4)+2 ,(dataBdim[2]*4)+2,1]== roww
#left 5,5,3 - testing right of one to the Left
roww = UInt32(0)
@setBitTo(roww,5,true)
@test  mainArr[(dataBdim[1]*4) ,(dataBdim[2]*4)+2,3]== roww
##right 5,5,4- testing left of block to the right
roww = UInt32(0)
@setBitTo(roww,5,true)
@test  mainArr[(dataBdim[1]*5)+1 ,(dataBdim[2]*4)+2,4]== roww

#posterior 5,5,5 - testring anterior of the  block to the back
roww = UInt32(0)
@setBitTo(roww,5,true)
@test mainArr[(dataBdim[1]*4)+2 ,(dataBdim[2]*4),5]== roww
#anterior 5,6,6 - testring posterior of the  block to the front
roww = UInt32(0)
@setBitTo(roww,5,true)
@test mainArr[(dataBdim[1]*4)+2 ,(dataBdim[2]*5)+1,6]== roww

########## checking results of those single voxels...

function checkIsInResList(resList,x,y,z,dir)::Bool
  println("looking for x $(x) y $(y) z $(z) dir $(dir)   [$(x),$(y),$(z),1,$(dir),0] ")
  for i in 1:size(resList)[1]
      if(resList[i,1]== x  && resList[i,2]== y && resList[i,3]== z && resList[i,5]== dir  )
        return true
      end
  end
  return false
end
cpuResList = Int64.(Array(resList))

cpuResList[261011,:]

filtered= []
for i in 1:size(resList)[1]
  if(cpuResList[i,1]>0)
    push!(filtered,Int64.(cpuResList[i,:]) )
  end
end
filtered

[130, 22, 32, 1, 5, 0]==[130,22,32,1,5,0]

#looking for [161,22,101,1,2,0] got

Int64.(Array(resList)[1,:])
@test checkIsInResList(cpuResList,(dataBdim[1]*4)+4 ,(dataBdim[2]*4)+4,(2-1)*dataBdim[3]+1, 6 )
@test checkIsInResList(cpuResList,(dataBdim[1]*4)+2 ,(dataBdim[2]*4)+2,(1-1)*dataBdim[3]+32, 5 )
@test checkIsInResList(cpuResList,(dataBdim[1]*4) ,(dataBdim[2]*4)+2,(3-1)*dataBdim[3]+5, 1 )

# @test checkIsInResList(cpuResList,(dataBdim[1]*5)+1 ,(dataBdim[2]*4)+2,(4-1)*dataBdim[3]+5, 2 )
@test checkIsInResList(cpuResList,(dataBdim[1]*4)+2 ,(dataBdim[2]*4),(5-1)*dataBdim[3]+5, 3)
@test checkIsInResList(cpuResList,(dataBdim[1]*4)+2 ,(dataBdim[2]*5)+1,(6-1)*dataBdim[3]+5, 4 )

# sum(resList)

afterDil = UInt32(0)
for bitPos in [1,2,4,5,6,31,32]
    @setBitTo(afterDil,bitPos,true)
end
afterDil
@test mainArr[33,dataBdim[2]+1,2]== afterDil
@test mainArr[33,dataBdim[2]*2,2]== afterDil
@test mainArr[64,dataBdim[2]+1,2]== afterDil

Int64(mainArr[33,dataBdim[2]+1,2])

rowOne = UInt32(0)
@setBitTo(rowOne,1,true)
@setBitTo(rowOne,5,true)
@setBitTo(rowOne,32,true)

  @test mainArr[33+1,dataBdim[2]+1,2]== rowOne
  @test mainArr[33+1,dataBdim[2]*2,2]== rowOne
  @test mainArr[64+1,dataBdim[2]+1,2]== rowOne

  @test mainArr[33-1,dataBdim[2]+1,2]== rowOne
  @test mainArr[33-1,dataBdim[2]*2,2]== rowOne
  @test mainArr[64-1,dataBdim[2]+1,2]== rowOne

for i in [-1,1]
  @test mainArr[33,dataBdim[2]+1+i,2]== rowOne
  @test mainArr[33,dataBdim[2]*2+i,2]== rowOne
  @test mainArr[64,dataBdim[2]+1+i,2]== rowOne
end
fullOnes = 0
for bitPos in 1:32
    @setBitTo(fullOnesTobecome,bitPos,true)
end

#this should be full after dilatation
for xx in ((3*32)+1):((4*32)), yy in ((3*dataBdim[2])+1):((4*dataBdim[2]))
  @test mainArr[xx,yy,3]== fullOnes
end

uu = Array(paddingStore)[2,2,2,:,:]
uu[1,1]
Int64(sum(uu))

#3) weather the information is block full has been properly set in the metadata the same is to be activated ...
 @test metaData[2,2,2,getFullInGoldNumb()]==0
 @test metaData[3,3,3,getFullInGoldNumb()]==1

 @test metaData[3,3,2,getIsToBeActivatedInGoldNumb()]==1
 @test metaData[3,3,4,getIsToBeActivatedInGoldNumb()]==1
 @test metaData[3,2,3,getIsToBeActivatedInGoldNumb()]==1
 @test metaData[3,4,3,getIsToBeActivatedInGoldNumb()]==1
 @test metaData[2,3,3,getIsToBeActivatedInGoldNumb()]==1
 @test metaData[4,3,3,getIsToBeActivatedInGoldNumb()]==1

#4)wheather in results  we have entries that should be present there so correct x,y,z and dir

@test checkIsInResList(cpuResList,(dataBdim[1])+1 ,(dataBdim[2])+1,dataBdim[3]+2,6)
@test checkIsInResList(cpuResList,(dataBdim[1])+1 ,(dataBdim[2])+1,dataBdim[3]+4,5)
@test checkIsInResList(cpuResList,(dataBdim[1]*2)-1 ,(dataBdim[2])+1, dataBdim[3]+5,1)
@test checkIsInResList(cpuResList,(dataBdim[1])+2 ,(dataBdim[2])+1,dataBdim[3]+5,2)
@test checkIsInResList(cpuResList,(dataBdim[1])*2 ,(dataBdim[2]*2)-1,dataBdim[3]+5,3)
@test checkIsInResList(cpuResList,(dataBdim[1])+1 ,(dataBdim[2])+2,dataBdim[3]+5,4)

# for xx in ((3*32)+1):((4*32)), yy in ((3*dataBdim[2])+1):((4*dataBdim[2])), zz in ((3*32)+1):((4*32))
#   if(iseven(zz))
#     #top is 6 and is checked before bottom
#     push!(shouldBeInResultSet,[xx,yy,zz,isGold,6,iterNumb] )
#   end
# end

# for entry in shouldBeInResultSet
#   @test checkIsInResList(resList,entry[1],entry[2],entry[3],entry[5])
# end

# using CUDA

# mean(CUDA.ones(8))

# # aa = mainArr[33,dataBdim[2]+1,2]
# aa= 16
# for i in 1:32
#   if(isBit1AtPos(aa,i))
#     print("i $(i)  ")
#   end
# end

# isBit1AtPos(2147483665,dataBdim[3])

# sum(paddingStore)
# if(xm ==2 && ym==2 && zm==2 && threadIdxX()==1 && threadIdxY()==1)     CUDA.@cuprint "7  locArr $(locArr) sourceShmem  $(shmemblockData[threadIdxX(),threadIdxY(),1]) res $(shmemblockData[threadIdxX(),threadIdxY(),2]) \n"     end

#5 check weather result counters are set to correct numbers

# #6) are the results in correct spots - weahter they are related to the ques the should be ...
# #where offset will be 0 for block 2,2,2 and 20000 for 3,3,3

# metaBlock = 0
# for xMetaa in 1:3,yMetaa in 1:3, zMetaa in 1:3
#   metaBlock+=20
#   for i in 1:14
#     metaData[2,2,2,getResOffsetsBeg()+metaBlock]=i*1000
#   end
# end

# function getBegEnd(dir)::Bool
#   for i in 1:14
#     if(i == )
#     beg = i*1000+twoTwoTwoOffset*15*1000
#     endd = (i+1)*1000+twoTwoTwoOffset*15*1000
#   end

# end

# function checkIsInResListAndSpot(resList,x,y,z,dir)::Bool
#   for i in 1:length(resList)
#       if(resList[i,:]== [x,y,z,isGold,dir,iterNumb])
#         return true
#       end
#   end
#   return false
# end

# for entry in shouldBeInResultSet
#   @test checkIsInResList(resList,entry[1],entry[2],entry[3],entry[5])
# end


text:
using AzManagers, Distributed, Documenter

makedocs(sitename="AzManagers", modules=[AzManagers])

deploydocs(
    repo = "github.com/ChevronETC/AzManagers.jl.git",
)

ext: jl


text:
"""
    function disagg!(d::Dict, set::Dict; save_build = true, overwrite = false)

# Arguments
- `d::Dict` of DataFrames containing the model data.
- `set::Dict` of Arrays describing region, sector, final demand, etc.

# Keywords
- `save_build = true`
- `overwrite = false`
See [`SLiDE.build`](@ref) for keyword argument descriptions.

# Returns
- `d::Dict` of DataFrames containing the model data at the disaggregation step
"""
function disagg(
    dataset::String,
    d::Dict,
    set::Dict;
    save_build::Bool = DEFAULT_SAVE_BUILD,
    overwrite::Bool = DEFAULT_OVERWRITE
    )
    CURR_STEP = "disagg"

    d = merge(d, Dict(
        :r => fill_with((r = set[:r],), 1.0),
        (:yr,:r,:g) => fill_with((yr = set[:yr], r = set[:r], g = set[:g]), 1.0)))

    # If there is already disaggregation data, read it and return.
    d_read = read_build(dataset, CURR_STEP; overwrite = overwrite)
    !(isempty(d_read)) && (return (d_read, set))

    # Run all disaggregation calculations.
    _set_gm!(d, set)

    d[:region] = edit_with(d[:region], Rename(:g,:s))
    _disagg_ys0!(d)
    _disagg_id0!(d)
    _disagg_ty0!(d, set)
    _disagg_va0!(d, set)
    _disagg_ld0!(d)
    _disagg_kd0!(d)

    d[:region] = edit_with(d[:region], Rename(:s,:g))
    _disagg_fdcat!(d)
    _disagg_g0!(d)
    _disagg_i0!(d)
    _disagg_cd0!(d)
    _disagg_c0!(d)

    d[:yh0_temp] = _disagg_yh0!(d)
    _disagg_fe0!(d)
    d[:x0_temp] = _disagg_x0!(d, set)
    d[:s0_temp] = _disagg_s0!(d)
    _disagg_a0!(d)

    _disagg_ta0!(d)
    _disagg_tm0!(d)

    _disagg_thetaa!(d)
    _disagg_m0!(d)
    _disagg_md0!(d)
    d[:rx0_temp] = _disagg_rx0!(d)

    _disagg_diff!(d)
    _apply_diff!(d, set)

    _disagg_bop!(d)
    _disagg_pt0!(d)
    _disagg_dc0!(d)

    _disagg_dd0!(d)
    _disagg_nd0!(d)

    _disagg_dm0!(d)
    _disagg_nm0!(d)
    _disagg_xd0!(d)
    _disagg_xn0!(d)
    _disagg_hhadj!(d)

    # Should other
    d[:xn0][d[:xn0][:,:value] .< 1e-8,:value] .= 0;
    d[:xd0][d[:xd0][:,:value] .< 1e-8,:value] .= 0;

    write_build!(dataset, CURR_STEP, d; save_build = save_build)
    write_build!(dataset, SET_DIR, Dict(k => set[k] for k in [:gm]))
    return (d, set)
end

"""
`ys(yr,r,s,g)`, regional sectoral output

```math
\\bar{ys}_{yr,r,s,g} = \\alpha_{yr,r,s}^{gsp} \\tilde{ys}_{yr,s,g}
```
"""
function _disagg_ys0!(d::Dict)
    println("  Disaggregating ys0(yr,r,s,g), regional sectoral output")
    :r in propertynames(d[:ys0]) && (return d[:ys0])

    d[:ys0] = d[:region] * d[:ys0]
    return d[:ys0]
end

"""
`id(yr,r,g,s)`, regional intermediate demand

```math
\\bar{id}_{yr,r,g,s} = \\alpha_{yr,r,s}^{gsp} \\tilde{id}_{yr,g,s}
```
"""
function _disagg_id0!(d::Dict)
    println("  Disaggregating id0(yr,r,g,s), regional intermediate demand")
    :r in propertynames(d[:id0]) && (return d[:id0])

    d[:id0] = d[:region] * d[:id0]
    return d[:id0]
end

"""
`ty(yr,r,s)`, production tax rate

```math
\\begin{aligned}
\\bar{ty}_{yr,r,s}^{rev} &= \\alpha_{yr,r,s}^{gsp} \\tilde{va}_{yr,va,s} \\;\\forall\\; va = othtax \\\\
\\bar{ty}_{yr,r,s} &= \\frac{\\tilde{ty}_{yr,r,s}}{\\sum_{g} \\bar{ys}_{yr,r,s,g}}
\\end{aligned}
```
"""
function _disagg_ty0!(d::Dict, set::Dict)
    println("  Disaggregating ty0(yr,r,s), production tax rate")
    :va in propertynames(d[:va0]) && _unstack_va0!(d, set)
    idx = findindex(d[:va0])

    ty0_rev = d[:region] * d[:va0][:,[idx;:othtax]]
    d[:ty0] = dropnan(ty0_rev / combine_over(d[:ys0], :g))
    return d[:ty0]
end

"""
`va(yr,va,s)`, regional value added

```math
\\bar{va}_{yr,r,s} = \\alpha_{yr,r,s}^{gsp} \\sum_{va = compen,surplus} \\tilde{va}_{yr,va,s}
```
"""
function _disagg_va0!(d::Dict, set::Dict)
    println("  Disaggregating va0(yr,r,s), regional share of value added.")
    :r in propertynames(d[:va0])  && (return d[:va0])
    :va in propertynames(d[:va0]) && (_unstack_va0!(d, set))
    idx = findindex(d[:va0])

    df = d[:va0][:,[idx;:compen]] + d[:va0][:,[idx;:surplus]]
    d[:va0] = d[:region] * df
    return d[:va0]
end

"""
`ld0(yr,r,s)`, labor demand

```math
\\bar{ld}_{yr,r,s} = \\theta_{yr,r,s}^{ls} \\bar{va}_{yr,s,g}
```
"""
function _disagg_ld0!(d::Dict)
    println("  Disaggregating ld0(yr,r,s), labor demand")
    !(:r in propertynames(d[:va0])) && _disagg_va0!(d, set)
    d[:ld0] = d[:labor] * d[:va0]
    return d[:ld0]
end

"""
`kd0(yr,r,s)`, capital demand

```math
\\bar{kd}_{yr,r,s} = \\bar{va}_{yr,r,s} - \\bar{ld}_{yr,r,s}
```
"""
function _disagg_kd0!(d::Dict)
    println("  Disaggregating kd0(yr,r,s), capital demand")
    d[:kd0] = d[:va0] - d[:ld0]
end

"""
    _disagg_fdcat!(d::Dict)
This function aggregates final demand categories into national consumption (`C`), government
(`G`), and investment (`I`) demand.
"""
function _disagg_fdcat!(d::Dict)
    if "pce" in d[:fd0][:,:fd]
        x = [
            Map(joinpath("crosswalk","fd.csv"), [:fd], [:fdcat], [:fd], [:fd], :inner),
            Combine("sum",propertynames(d[:fd0])),
        ]
        d[:fd0] = edit_with(d[:fd0], x)
        # d[:fd0] = combine_over(d[:fd0], :fd)
    end
    return d[:fd0]
end

"""
`g0(yr,r,g)`, national government demand

```math
\\bar{g}_{yr,r,g} = \\alpha_{yr,r,g}^{sgf} \\sum_{G \\in fd} \\tilde{fd}_{yr,g,fd}
```
"""
function _disagg_g0!(d::Dict)
    println("  Disaggregating g0(yr,r,g), national government demand")
    ("pce" in d[:fd0][:,:fd]) && (_disagg_fdcat!(d))

    df = filter_with(d[:fd0], (fd="G",); drop = true)
    d[:g0] = d[:sgf] * df
    return d[:g0]
end

"""
`i0(yr,r,g)`, national investment demand

```math
\\bar{i}_{yr,r,g} = \\alpha_{yr,r,g}^{gsp} \\sum_{I \\in fd} \\tilde{fd}_{yr,g,fd}
```
"""
function _disagg_i0!(d::Dict)
    println("  Disaggregating i0(yr,r,g), national investment demand")
    !("pce" in d[:fd0][:,:fd]) && (_disagg_fdcat!(d))

    df = filter_with(d[:fd0], (fd="I",); drop = true)
    d[:i0] = d[:region] * df
    return d[:i0]
end

"""
`cd0(yr,r,g)`, national final consumption

```math
\\bar{cd}_{yr,r,g} = \\alpha_{yr,r,g}^{pce} \\sum_{C \\in fd} \\tilde{fd}_{yr,g,fd}
```
"""
function _disagg_cd0!(d::Dict)
    println("  Disaggregating cd0(yr,r,g), national final consumption")
    !("pce" in d[:fd0][:,:fd]) && (_disagg_fdcat!(d))

    df = filter_with(d[:fd0], (fd="C",); drop = true)
    d[:cd0] = d[:pce] * df
    return d[:cd0]
end

"""
`c0(yr,r)`, total final household consumption

```math
\\bar{c}_{yr,r} = \\sum_{g} \\bar{cd}_{yr,r,g}
```
"""
function _disagg_c0!(d::Dict)
    println("  Disaggregating c0(yr,r), total final household consumption")
    !(:cd0 in keys(d)) && _disagg_cd0!(d)
    d[:c0] = combine_over(d[:cd0], :g)
    return d[:c0]
end

"""
`yh0(yr,r,g)`, household production

```math
\\bar{yh}_{yr,r,g} = \\alpha_{yr,r,g} \\tilde{fs}_{yr,g}
```
"""
function _disagg_yh0!(d::Dict)
    if !(:diff in keys(d))
        println("  Disaggregating yh0(yr,r,g), household production")
        d[:yh0] = d[:region] * d[:fs0]
    else
        println("  Applying difference to yh0(yr,r,g), household production")
        d[:yh0] = d[:yh0] + d[:diff]
    end
    return dropmissing!(d[:yh0])
end

"""
`fe0(yr,r)`, total factor supply

```math
\\bar{fe}_{yr,r} = \\sum_{s} \\bar{va}_{yr,r,s}
```
"""
function _disagg_fe0!(d::Dict)
    println("  Disaggregating fe0(yr,r), total factor supply")
    d[:fe0] = combine_over(d[:va0], :s)
    return d[:fe0]
end

"""
`x0(yr,r,g)`, foreign exports

```math
\\bar{x}_{yr,r,g} = \\alpha_{yr,r,g}^{utd} \\tilde{x}_{yr,g}
```
"""
function _disagg_x0!(d::Dict, set::Dict)
    if !(:diff in keys(d))
        println("  Disaggregating x0(yr,r,g), foreign exports")
        !(:notrd in keys(set)) && _set_notrd!(d, set)

        df_exports = filter_with(d[:utd], (t = "exports",); drop = true)
        df_region = filter_with(d[:region], (g = set[:notrd],))

        df_trd = dropmissing(df_exports * d[:x0])
        df_notrd = dropmissing(df_region * d[:x0])

        d[:x0] = [df_trd; df_notrd]
    else
        println("  Applying difference to x0(yr,r,g), foreign exports")
        d[:x0] = d[:x0] + d[:diff]
    end

    return dropmissing!(d[:x0])
end

"""
`s0(yr,r,g)`, total supply

```math
\\bar{s}_{yr,r,g} = \\sum_{s} \\bar{ys}_{yr,r,s,g} + \\bar{yh}_{yr,r,g}
```
"""
function _disagg_s0!(d::Dict)
    if !(:diff in keys(d))
        println("  Disaggregating s0(yr,r,g), total supply")
        d[:s0] = combine_over(d[:ys0], :s) + d[:yh0]
    else
        println("  Applying difference to s0(yr,r,g), total supply")
        d[:s0] = d[:s0] + d[:diff]
    end
    return dropmissing!(d[:s0])
end

"""
`a0(yr,r,g)`, domestic absorption

```math
a_{yr,r,g} = \\bar{cd}_{yr,r,g} + \\bar{g}_{yr,r,g} + \\bar{i}_{yr,r,g} + \\sum_{s}\\bar{id}_{yr,r,g}
```
"""
function _disagg_a0!(d::Dict)
    println("  Disaggregating a0(yr,r,g), domestic absorption")

    d[:a0] = dropmissing(d[:cd0] + d[:g0] + d[:i0] + combine_over(d[:id0], :s))
    return d[:a0]
end

"`ta0(yr,r,g)`: Absorption taxes"
function _disagg_ta0!(d::Dict)
    println("  Disaggregating ta0(yr,r,g), absorption taxes")
    d[:ta0] = d[:ta0] * d[:r]
    return d[:ta0]
end

"`tm0(yr,r,g)`: Import taxes"
function _disagg_tm0!(d::Dict)
    println("  Disaggregating tm0(yr,r,g), import taxes")
    d[:tm0] = d[:tm0] * d[:r]
    return d[:tm0]
end

"""
`thetaa(yr,r,g)`, share of regional absorption

```math
\\alpha_{yr,r,g}^{abs} = \\frac{\\bar{a}_{yr,r,g}}{\\sum_{rr}\\bar{a}_{yr,r,g}}
```
"""
function _disagg_thetaa!(d::Dict)
    println("  Disaggregating thetaa(yr,r,g), share of regional absorption")
    d[:thetaa] = dropnan(d[:a0] / transform_over(d[:a0], :r))
    return d[:thetaa]
end

"""
`m0(yr,r,g)`, foreign Imports

```math
\\bar{m}_{yr,r,g} = \\alpha_{yr,r,g}^{abs} \\tilde{m}_{yr,g}
```
"""
function _disagg_m0!(d::Dict)
    println("  Disaggregating m0(yr,r,g), foreign imports")
    d[:m0] = d[:thetaa] * d[:m0]
    return d[:m0]
end

"""
`md0(yr,r,m,g)`, margin demand

```math
\\bar{md}_{yr,r,m,g} = \\alpha_{yr,r,g}^{abs} \\tilde{md}_{yr,m,g}
```
"""
function _disagg_md0!(d::Dict)
    println("  Disaggregating md0(yr,r,m,g), margin demand")
    d[:md0] = dropmissing(d[:thetaa] * d[:md0])
    return d[:md0]
end

"""
`rx0(yr,r,g)`, re-exports

```math
\\bar{rx}_{yr,r,g} = \\bar{x}_{yr,r,g} - \\bar{s}_{yr,r,g}
```
"""
function _disagg_rx0!(d::Dict; round_digits = DEFAULT_ROUND_DIGITS)
    if !(:diff in keys(d))
        println("  Disaggregating rx0(yr,r,g), re-exports")
        d[:rx0] = d[:x0] - d[:s0]

        if round_digits !== false
            d[:rx0][round.(d[:rx0][:,:value]; digits = round_digits) .< 0, :value] .= 0.0
        end
    else
        println("  Applying difference to rx0(yr,r,g), re-exports")
        d[:rx0] = d[:rx0] + d[:diff]
    end

    return dropmissing(d[:rx0])
end

"""
`dc0`,

```math
\\bar{dc}_{yr,r,g} = \\bar{s}_{yr,r,g} - \\bar{x}_{yr,r,g} + \\bar{rx}_{yr,r,g}
```
"""
function _disagg_dc0!(d::Dict; round_digits = DEFAULT_ROUND_DIGITS)
    # (!!!!) name for this?
    d[:dc0] = dropmissing((d[:s0] - d[:x0] + d[:rx0]))

    # if round_digits !== false
    #     d[:dc0][!,:value] .= round.(d[:dc0][:,:value]; digits = round_digits)
    # end
    (round_digits !== false) && (d[:dc0] = round!(d[:dc0], :value; digits = round_digits))
    return d[:dc0]
end

"""
`pt0`,

```math
\\begin{aligned}
\\bar{pt}_{yr,r,g} = &\\left(1 - \\bar{ta}_{yr,r,g} \\right) \\bar{a}_{yr,r,g} + \\bar{rx}_{yr,r,g}
\\\\               - &\\left(1 + \\bar{tm}_{yr,r,g} \\right) \\bar{m}_{yr,r,g} - \\sum_{m} \\bar{md}_{yr,r,m,g}
\\end{aligned}
```
"""
function _disagg_pt0(d::Dict; round_digits = DEFAULT_ROUND_DIGITS)
    df_pta0 = dropmissing(((d[:yr,:r,:g] - d[:ta0]) * d[:a0]) + d[:rx0])
    df_ptm0 = dropmissing(((d[:yr,:r,:g] + d[:tm0]) * d[:m0]) + combine_over(d[:md0], :m))

    df = df_pta0 - df_ptm0
    (round_digits !== false) && round!(df, :value; digits = round_digits)
    return df
end

function _disagg_pt0!(d::Dict; round_digits = DEFAULT_ROUND_DIGITS)
    d[:pt0] = _disagg_pt0(d, round_digits = round_digits)
    return d[:pt0]
end

"`diff`:"
function _disagg_diff!(d::Dict; round_digits = DEFAULT_ROUND_DIGITS)
    # df = _disagg_pt0!(d)
    df = _disagg_pt0(d; round_digits = DEFAULT_ROUND_DIGITS)

    df[!,:value] .= - min.(0, df[:,:value])
    d[:diff] = dropzero!(df)
end

"""
Add an adjustment to `rx_{yr,r,g}`, `s_{yr,r,g}`, `x_{yr,r,g}`, and `yh_{yr,r,g}`.
"""
function _apply_diff!(d::Dict, set::Dict)
    _disagg_rx0!(d)
    _disagg_s0!(d)
    _disagg_x0!(d, set)
    _disagg_yh0!(d)
end

"""
`bopdef0(yr,r)`, balance of payments (closure parameter)

```math
\\bar{bop}_{yr,r} = \\sum_{g} \\left( \\bar{m}_{yr,r,g} - \\bar{x}_{yr,r,g} \\right)
```
"""
function _disagg_bop!(d::Dict)
    println("  Disaggregating bopdef0(yr,r), balance of payments (closure parameter)")
    d[:bopdef0] = combine_over((d[:m0] - d[:x0]), :g)
    return d[:bopdef0]
end

"`gm`: Commodities employed in margin supply"
function _set_gm!(d::Dict, set::Dict)
    ms0_sum = combine_over(d[:ms0], [:yr,:m])
    md0_sum = combine_over(d[:md0], [:yr,:m])

    ms0_sum[!,:value] .= ms0_sum[:,:value] .!= 0.0
    md0_sum[!,:value] .= md0_sum[:,:value] .!= 0.0

    gm = ms0_sum + md0_sum
    set[:gm] = gm[gm[:,:value] .> 0, :g]
    return set[:gm]
end

"""
`dd0max(yr,r,g)`, maximum regional demand from local market

```math
\\hat{dd}_{yr,r,g} = \\min\\left\\{\\bar{pt}_{yr,r,g}, \\bar{dc}_{yr,r,g} \\right\\}
```
"""
function _disagg_dd0max!(d::Dict)
    println("  Disaggregating dd0max(yr,r,g), maximum regional demand from local market")
    cols = propertynames(d[:pt0])
    df = indexjoin(d[:pt0], d[:dc0]; id = [:pt0,:dc0])
    df[!,:value] = min.(df[:,:pt0], df[:,:dc0])

    d[:dd0max] = df[:,cols]
    return d[:dd0max]
end

"`nd0max(yr,r,g)`, maximum regional demand from national market"
function _disagg_nd0max!(d::Dict)
    println("  Disaggregating nd0max(yr,r,g), maximum regional demand from national market")
    cols = propertynames(d[:pt0])
    df = indexjoin(d[:pt0], d[:dc0]; id = [:pt0,:dc0])
    df[!,:value] = min.(df[:,:pt0], df[:,:dc0])

    d[:nd0max] = df[:,cols]
    return d[:nd0max]
end

"`dd0min(yr,r,g)`, minimum regional demand from local market"
_disagg_dd0min(d::Dict) = d[:pt0] - d[:dd0max]

"`nd0min(yr,r,g)`, minimum regional demand from national market"
_disagg_nd0min(d::Dict) = d[:pt0] - d[:nd0max]

"""
`dd0(yr,r,g)`, regional demand from local market

```math
\\bar{dd}_{yr,r,g} = \\rho_{r,g}^{cfs} \\hat{dd}_{yr,r,g}
```
"""
function _disagg_dd0!(d::Dict)
    println("  Disaggregating dd(yr,r,g), regional demand from local market")
    !(:dd0max in keys(d)) && _disagg_dd0max!(d)
    d[:dd0] = d[:dd0max] * d[:rpc]
    return d[:dd0]
end

"""
`nd0_(yr,r,g)`, regional demand from national market

```math
\\bar{nd}_{yr,r,g} = \\bar{pt}_{yr,r,g} - \\bar{dd}_{yr,r,g}
```
"""
function _disagg_nd0!(d::Dict; round_digits = DEFAULT_ROUND_DIGITS)
    println("  Disaggregating nd(yr,r,g), regional demand from national market")
    df_pt0 = _disagg_pt0(d; round_digits = false)

    d[:nd0] = df_pt0 - d[:dd0]
    # (round_digits !== false) && (d[:nd0][!,:value] .= round.(d[:nd0][:,:value]; digits = round_digits))
    (round_digits !== false) && (d[:nd0] = round!(d[:nd0], :value; digits = round_digits))
    return d[:nd0]
end

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

"""
`mrgshr(yr,r,m)`, share of margin demand by region

```math
\\alpha_{yr,r,m}^{md} = \\frac{\\sum_{g}\\bar{md}_{yr,r,m,g}}{\\sum_{r,g}\\bar{md}_{yr,r,m,g}}
```
"""
function _disagg_mrgshr(d::Dict)
    # (!!!!) notation alpha^{md} to match thetaa, alpha^a
    df = combine_over(d[:md0], :g)
    df = df / transform_over(df, :r)
    return df
end

"""
`ms0tot(yr,r,m,g)`, designate total supply of margins

```math
\\hat{ms}_{yr,r,m,g} = \\alpha_{yr,r,m}^{md} \\bar{ms}_{yr,g,m}
```
"""
function _disagg_ms0tot!(d::Dict)
    d[:ms0tot] = _disagg_mrgshr(d) * d[:ms0]
    return d[:ms0tot]
end

"""
`shrtrd(yr,r,g,m)`, share of margin total by margin type

```math
\\beta_{yr,r,g,m}^{mar} = \\frac{\\hat{ms}_{yr,r,g,m}}{\\sum_{m}\\hat{ms}_{yr,r,g,m}}
```
"""
function _disagg_shrtrd!(d::Dict)
    !(:ms0tot in keys(d)) && _disagg_ms0tot!(d)

    df = d[:ms0tot]
    d[:shrtrd] = dropnan(df / transform_over(df, :m))
    return d[:shrtrd]
end

"""
`dm0(yr,r,g,m)`, margin supply from the local market

```math
\\bar{dm}_{yr,r,g,m} = \\min\\left\\{ \\rho_{r,g}^{cfs}\\hat{ms}_{yr,r,g,m},
    \\beta_{yr,r,m,g}^{mar} \\left(\\bar{dc}_{yr,r,g} - \\bar{dd}_{yr,r,g}\\right) \\right\\}
```
"""
function _disagg_dm0!(d::Dict)
    println("  Disaggregating dm0(yr,r,g,m), margin supply from the local market")
    !(:ms0tot in keys(d)) && _disagg_ms0tot!(d)
    !(:shrtrd in keys(d)) && _disagg_shrtrd!(d)

    cols = propertynames(d[:ms0tot])
    dm1 = dropmissing(d[:ms0tot] * d[:rpc])
    dm2 = dropmissing((d[:shrtrd] * d[:dc0]) - d[:dd0])

    df = indexjoin(dm1, dm2; id = [:dm1,:dm2])
    df[!,:value] .= min.(df[:,:dm1], df[:,:dm2])
    d[:dm0] = df[:,cols]
    return d[:dm0]
end

"""
`nm0(yr,r,g,m)`, margin demand from the national market

```math
\\bar{nm}_{yr,r,g,m} = \\hat{ms}_{yr,r,g,m} - \\bar{dm}_{yr,r,g,m}
```
"""
function _disagg_nm0!(d::Dict)
    println("  Disaggregating nm0(yr,r,g,m), margin demand from the national market")
    d[:nm0] = d[:ms0tot] - d[:dm0]
    return d[:nm0]
end

"""
`xd0(yr,r,g)`, regional supply to local market

```math
\\bar{xd}_{yr,r,g} = \\sum_{m}\\bar{dm}_{yr,r,g,m} + \\bar{dd}_{yr,r,g}
```
"""
function _disagg_xd0!(d::Dict)
    println("  Disaggregating xd0(yr,r,g), regional supply to local market")
    d[:xd0] = combine_over(d[:dm0], :m) + d[:dd0]
    return d[:xd0]
end

"""
`xn0(yr,r,g)`, regional supply to national market

```math
\\bar{xn}_{yr,r,g} = \\bar{s}_{yr,r,g} + \\bar{rx}_{yr,r,g} - \\bar{xd}_{yr,r,g} - \\bar{x}_{yr,r,g}
```
"""
function _disagg_xn0!(d::Dict)
    println("  Disaggregating xn0(yr,r,g), regional supply to national market")
    d[:xn0] = d[:s0] + d[:rx0] - d[:xd0] - d[:x0]
    return d[:xn0]
end

"""
`hhadj(yr,r)`, household adjustment

```math
\\begin{aligned}
\\bar{adj^{hh}}_{yr,r} = &\\bar{c}_{yr,r}
\\\\ &- \\sum_{s}\\left( \\bar{ld}_{yr,r,s} + \\bar{kd}_{yr,r,s} + \\bar{yh}_{yr,r,s} \\right) - \\bar{bop}_{yr,r}
\\\\ &- \\sum_{s}\\left( \\bar{ta}_{yr,r,s}\\bar{a0}_{yr,r,s} + \\bar{tm}_{yr,r,s}\\bar{m}_{yr,r,s} + \\bar{ty}_{yr,r,s}\\sum_{g}\\bar{ys}_{yr,r,s,g} \\right)
\\\\ &+ \\sum_{s}\\left( \\bar{g}_{yr,r,s} + \\bar{i}_{yr,r,s} \\right)
\\end{aligned}
```
"""
function _disagg_hhadj!(d::Dict)
    println("  Disaggregating hhadj(yr,r), household adjustment")
    dh = Dict(k => edit_with(copy(d[k]), Rename(:g,:s))
        for k in [:c0,:ld0,:kd0,:yh0,:bopdef0,:ta0,:a0,:tm0,:ty0,:g0,:i0,:m0])
    dh[:ys0] = copy(d[:ys0])

    d[:hhadj] = dh[:c0] -
        combine_over(dh[:ld0] + dh[:kd0] + dh[:yh0], :s) - dh[:bopdef0] -
        combine_over(dh[:ta0]*dh[:a0] + dh[:tm0]*dh[:m0] + dh[:ty0]*combine_over(dh[:ys0],:g), :s) +
        combine_over(dh[:g0] + dh[:i0], :s)

    return d[:hhadj]
end

function _unstack_va0!(d::Dict, set::Dict)
    if :va in propertynames(d[:va0])
        d[:va0] = edit_with(unstack(d[:va0], :va, :value),
            Replace.(Symbol.(set[:va]), missing, 0.0))
    end
    return d[:va0]
end

max_stars_repo_name: NREL/SLiDE

max_issues_repo_name: NREL/SLiDE

max_forks_repo_name: NREL/SLiDE


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "lrslib"
version = v"0.2.0"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/JuliaPolyhedra/lrslib.git",
              "d8b723a2c315614979a8354f9e768d273d14a215"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/lrslib*
extraargs=""
cflags="-fPIC -O3 -Wall"

# 32bit linux, arm and windows:
if [[ $target == i686* ]] || [[ $target == arm* ]]; then
  # no 128 bit ... patch makefile
  atomic_patch -p1 ../patches/no128bit.patch
fi

if [[ $target == *apple* ]]; then
  export CC=gcc
  sed -i -e 's#-Wl,-soname=#-install_name #' makefile
  extraargs="SONAME=liblrs.0.dylib SHLINK=liblrs.dylib SHLIB=liblrs.0.0.0.dylib"
elif [[ $target == *freebsd* ]]; then
  export CC="gcc"
elif [[ $target == *mingw* ]]; then
  extraargs="SONAME=liblrs-0.dll SHLINK=liblrs.dll SHLIB=liblrs-0-0-0.dll"
  cflags="$cflags -DSIGNALS -DTIMES"
fi

make prefix=$prefix INCLUDEDIR=$prefix/include LIBDIR=${libdir} CFLAGS="$cflags" $extraargs -j ${nproc} install
if [[ $target == *mingw* ]]; then
  # rename binaries and move libraries
  for file in ${bindir}/{lrs,lrsnash,redund}; do mv $file $file.exe; done
  mv ${prefix}/lib/*lrs*.dll ${libdir}/
fi

${CC} -shared ${cflags} -o "${libdir}/liblrsnash.${dlext}" lrsnashlib.c -L${libdir} -llrs -lgmp -Wl,-rpath,${libdir} -DMA -DGMP -DLRS_QUIET -I${includedir}
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    ExecutableProduct("lrs", :lrs)
    ExecutableProduct("lrsnash", :lrsnash)
    ExecutableProduct("redund", :redund)
    LibraryProduct("liblrs", :liblrs)
    LibraryProduct("liblrsnash", :liblrsnash)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("GMP_jll", v"6.1.2"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: waralex/Yggdrasil

max_issues_repo_name: waralex/Yggdrasil

max_forks_repo_name: waralex/Yggdrasil


text:
"""
    norm_data(A; dims=1, use_MAD::Bool=false, corrected::Bool=false)

Computes ``\\frac{A - \\bar A}{\\sigma}``, where `dims` is 1 (the default) if `A` is
column-oriented and 2 if it is row-oriented. If `use_MAD` is `true`, ``\\bar A`` is
`median(A, dims=dims)` and ``\\sigma`` is the MAD version of the standard deviation. If
`use_MAD` is `false`, ``\\bar A`` is `mean(A, dims=dims)` and ``\\sigma`` is the standard
deviation. If `corrected` is `true` then the sum is scaled with `n-1`, whereas the sum is
scaled with `n` if `corrected` is `false`, where `n = size(A, dims)`.
"""
function norm_data(A; dims=1, use_MAD=false, corrected=false)
    A = Matrix{Float64}(A)
    if use_MAD
        return (A .- median(A, dims=dims)) ./ mapslices(mad, A, dims=dims)
    else
        return (A .- mean(A, dims=dims)) ./ std(A, dims=dims, corrected=corrected)
    end
end

"""
    classic_PCA(A; dims=1, standardize::Bool=true, corrected::Bool=false)

Performs PCA on a copy of the ``m × n`` matrix `A`. `dims` is 1 (the default) if `A` is
column-oriented and 2 if it is row-oriented. If `standardize` is `true` (the default), the
copy will be transformed using the MAD version of the standard deviation. If `corrected`
is `true` then the sum is scaled with `n-1`, whereas the sum is scaled with `n` if
`corrected` is `false`, where `n = size(A, dims)`.

Returns `(V, S, prod)`, where ``A = U Σ V^*``, with ``Σ = \\text{diagm}(S)`` being the
diagonal matrix of singular values sorted in descending order and `prod` is the matrix
multiplication of `V` and `A`, with `A` being standardized by default.
"""
function classic_PCA(A; dims=1, standardize=true, corrected=false)
    data = copy(A)
    data = standardize ? norm_data(data, dims=dims, use_MAD=true, corrected=corrected) : data
    U, S, V = svd(cov(data, dims=dims, corrected=corrected))
    prod = dims == 1 ? data * V : V' * data
    return V, S, prod
end


text:
module DynamicForestBehaviors

using AutomotiveDrivingModels

using RandomForestBehaviors.MvDecisionTrees
using Distributions.MvNormal

import AutomotiveDrivingModels:
    AbstractVehicleBehavior,
    AbstractVehicleBehaviorTrainParams,
    AbstractVehicleBehaviorPreallocatedData,

    preallocate_learning_data,
    select_action,
    calc_action_loglikelihood,
    train

export
    DynamicForestBehavior,
    DF_TrainParams,
    AutoregressiveMvNormLeaf,

    select_action,
    calc_action_loglikelihood,
    train

global GLOBAL_DIAG_SHRINKAGE_THRESOLD = 1e-6
DEFAULT_AUTOREGRESSION_CONSTANT = 0.5
DEFAULT_NUM_AUTOREGRESSION_PREDICTORS = 3
DEFAULT_NUM_PREDICTOR_SAMPLES = 10

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

immutable AutoregressiveMvNormLeaf <: MvDecisionTrees.Leaf
    A::Matrix{Float64} # regression matrix A = UΦᵀ(ΦΦᵀ + γI)⁻¹ where γ is a ridge regularization param, U is column-wise targets, Φ is column-wise features
    m::MvNormal        # μ = Aϕ must be overwritten on each use
    predictor_indeces::Vector{Int} # indeces of predictors
end

function _regress_on_predictor_subset{T<:AbstractFloat, S<:Real}(
    data::MvDecisionTrees.TreeData{T,S},
    assignment_id::Int,
    Γ::Matrix{T}, # ϕ+1 × ϕ+1 diagonal γ matrix (preallocated)
    U::Matrix{T}, #   o × m
    Φ::Matrix{T}, # ϕ+1 × m (is overwritten)
    A::Matrix{T}, #   o × ϕ+1 (is overwritten)
    predictor_indeces::Vector{Int}, # length ϕ
    )

    #=
    Returns the loss when using this predictor set
    Overwrites Φ
    =#

    n = size(data.X,1)
    o, m = size(U)
    ϕ = length(predictor_indeces)
    @assert(size(Φ,1) == ϕ+1)
    @assert(size(Φ,2) == m)
    @assert(size(Γ,1) == size(Γ,2) == ϕ+1 )

    u_index = 0
    for i in 1 : n
        if data.assignment[i] == assignment_id
            u_index += 1
            for (j,k) in enumerate(predictor_indeces)
                Φ[j,u_index] = data.X[i,k] # features : [nrows, nfeatures], Φ: [nfeatures, nrows]
            end
            Φ[ϕ+1,u_index] = 1.0 # constant 1-predictor
        end
    end

    # regress

    A[:] = U*Φ'/(Φ*Φ' + Γ)

    # calc residual loss (square deviation)

    loss = 0.0
    for i in 1 : n
        if data.assignment[i] == assignment_id
            for j in 1 : o
                r = data.Y[j,i] - A[j,ϕ+1]
                for k in 1 : ϕ
                    r -= A[j,k]*data.X[i,predictor_indeces[k]]
                end
                loss += r*r # squared error
            end
        end
    end

    loss
end
function MvDecisionTrees.build_leaf{T<:AbstractFloat, S<:Real}(::Type{AutoregressiveMvNormLeaf},
    data::MvDecisionTrees.TreeData{T,S},
    assignment_id::Int,
    γ::T=DEFAULT_AUTOREGRESSION_CONSTANT, # ridge-regularization parameter
    n_autoregression_predictors::Int = DEFAULT_NUM_AUTOREGRESSION_PREDICTORS,
    n_random_predictor_samples::Int = DEFAULT_NUM_PREDICTOR_SAMPLES,
    )

    #=
    Fit the regression matrix A = UΦᵀ(ΦΦᵀ + γI)⁻¹
       γ is a ridge regularization param ≥ 0
         0 → least squares
         large γ → increased penalty on regression coefficients
       U is column-wise targets
       Φ is column-wise features plus a constant 1-feature

    This function picks a set of 'n_autoregression_predictors' predictors used to regress on the targets
    This is done 'n_random_predictor_samples' times
    The best in terms of residual error is kept
    The covariance matrix is static and computed based on the final residual errors
    =#

    @assert(n_autoregression_predictors ≥ 0)
    @assert(n_random_predictor_samples > 0)

    n, o, p = size(data)
    id_count = MvDecisionTrees.get_id_count(data, assignment_id)
    @assert(id_count > 0)

    ##################################
    # pull target matrix

    U = Array(T, o, id_count) # column-wise contatenation of target variables

    u_index = 0
    for i in 1 : n
        if data.assignment[i] == assignment_id
            u_index += 1
            for j in 1 : o
                U[j,u_index] = data.Y[j,i]
            end
        end
    end

    #################################
    # pull predictor matrix

    ϕ = min(n_autoregression_predictors, p) # number of real features
    Φ = Array(T, ϕ+1, id_count) # column-wise contatenation of predictors (also including constant 1-predictor)
    Γ = diagm(fill(γ, (ϕ+1)))
    A = Array(T, o, ϕ+1)
    predictor_indeces = Array(Int, ϕ)

    if ϕ == 0
        _regress_on_predictor_subset(data, assignment_id, Γ, U, Φ, A, predictor_indeces) # just solve for A
    else
        # pick several samples of predictor_indeces
        temp_A = deepcopy(A)
        temp_predictor_indeces = deepcopy(predictor_indeces)
        best_loss = Inf

        for i = 1 : n_random_predictor_samples
            MvDecisionTrees._reservoir_sample!(temp_predictor_indeces, p, ϕ)
            loss = _regress_on_predictor_subset(data, assignment_id, Γ, U, Φ, temp_A, temp_predictor_indeces)
            if loss < best_loss
                best_loss = loss
                copy!(A, temp_A)
                copy!(predictor_indeces, temp_predictor_indeces)
            end
        end

        # in case we found nothing
        if isinf(best_loss)
            fill!(A, 0.0)
            fill!(predictor_indeces, 1)
        end
    end

    ###########################################################3
    # calc covariance

    Σ = zeros(T, o, o)
    r_vec = Array(Float64, o)

    # obtain upper triangle
    for i in 1 : n
        if data.assignment[i] == assignment_id

            for j in 1 : o
                r_vec[j] = data.Y[j,i] - A[j,ϕ+1]
                for k in 1 : ϕ
                    r_vec[j] -= A[j,k]*data.X[i,predictor_indeces[k]]
                end
            end

            for a in 1 : o
                for b in a : o
                    Σ[a,b] += r_vec[a] * r_vec[b]
                end
            end
        end
    end

    if id_count < 2
        warn("Not enough samples to build a covariance matrix")
        for i in 1 : o
            Σ[i,i] = 1.0
        end
    else
        for i in 1 : o*o
            @assert(!isnan(Σ[i]))
            Σ[i] /= (id_count-1)
        end

        ε = GLOBAL_DIAG_SHRINKAGE_THRESOLD
        MvDecisionTrees._diagonal_shrinkage!(Σ, ε)
        for a in 2:o
            for b in 1:a-1
                Σ[a,b] = Σ[b,a]
            end
        end
    end

    AutoregressiveMvNormLeaf(A, MvNormal(Array(T, o), Σ), predictor_indeces)
end

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

type DF_TrainParams <: AbstractVehicleBehaviorTrainParams

    targets::ModelTargets
    indicators::Vector{AbstractFeature}

    ntrees::Int
    max_tree_depth::Int
    n_split_tries::Int
    min_samples_split::Int
    min_samples_leaves::Int
    n_random_predictor_samples::Int
    n_autoregression_predictors::Int
    use_PCA::Bool

    diag_shrinkage_threshold::Float64
    min_split_improvement::Float64
    partial_sampling::Float64
    autogression_coef::Float64

    function DF_TrainParams(;
        targets::ModelTargets = ModelTargets(Features.FUTUREDESIREDANGLE, Features.FUTUREACCELERATION),
        indicators::Vector{AbstractFeature} = AbstractFeature[],

        ntrees::Integer=10,
        max_tree_depth::Integer=5,
        n_split_tries::Integer=50,
        min_samples_split::Integer=100,
        min_samples_leaves::Integer=20,
        n_random_predictor_samples::Integer=10,
        n_autoregression_predictors::Integer=2,
        use_PCA::Bool=false,

        diag_shrinkage_threshold::Float64=1e-6,
        min_split_improvement::Float64=0.0,
        partial_sampling::Float64=0.7,
        autogression_coef::Float64=0.1,
        )

        retval = new()
        retval.targets = targets
        retval.indicators = indicators
        retval.ntrees = ntrees
        retval.max_tree_depth = max_tree_depth
        retval.n_split_tries = n_split_tries
        retval.min_samples_split = min_samples_split
        retval.min_samples_leaves = min_samples_leaves
        retval.n_random_predictor_samples = n_random_predictor_samples
        retval.n_autoregression_predictors = n_autoregression_predictors
        retval.use_PCA = use_PCA

        retval.diag_shrinkage_threshold = diag_shrinkage_threshold
        retval.min_split_improvement = min_split_improvement
        retval.partial_sampling = partial_sampling
        retval.autogression_coef = autogression_coef
        retval
    end
end
type DF_PreallocatedData <: AbstractVehicleBehaviorPreallocatedData

    # the full dataset matrices pre-extracted
    X::Matrix{Float64} # [nfeatures × nframes]
    Y::Matrix{Float64} # [2 × nframes]
    Z::Matrix{Float64} # [n_pca_features × nframes]

    extractor::FeatureSubsetExtractor
    preprocess::DataPreprocessor
    preprocess_pca::DataPreprocessor
    action_clamper::DataClamper

    function DF_PreallocatedData(dset::ModelTrainingData2, params::DF_TrainParams)

        retval = new()

        targets = params.targets
        indicators = params.indicators
        trainingframes = dset.dataframe
        nframes = nrow(trainingframes)
        nindicators = length(indicators)

        X = Array(Float64, nindicators, nframes) # NOTE(tim): initially indicators × nframes, we later transpose it
        Y = Array(Float64, 2, nframes)

        pull_design_and_target_matrices!(X, Y, trainingframes, targets, indicators)

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

        extractor = FeatureSubsetExtractor(indicators)
        preprocess = ChainedDataProcessor(extractor)
        push!(preprocess, X, DataNaReplacer, extractor.indicators)
        process!(X, preprocess.processors[end]) # process in place
        push!(preprocess, X, DataScaler) # ensure data is standardized before running PCA
        process!(X, preprocess.processors[end]) # process in place
        retval.X = copy(X)

        # Add a PCA scaler
        preprocess_pca = deepcopy(preprocess)
        n_PCA_features = min(nindicators, 20)
        push!(preprocess_pca, X, DataLinearTransform, n_PCA_features)

        Z = Array(Float64, n_PCA_features, nframes)
        retval.Z = process!(Z, X, preprocess_pca.processors[end])

        retval.Y = Y
        retval.extractor = extractor
        retval.preprocess = preprocess
        retval.preprocess_pca = preprocess_pca
        retval.action_clamper = DataClamper(
                Array(Float64, 2),
                vec(minimum(Y, 2)),
                vec(maximum(Y, 2))
            )

        retval
    end
end
type DynamicForestBehavior <: AbstractVehicleBehavior
    forest::Ensemble

    targets::ModelTargets
    indicators::Vector{AbstractVehicleBehavior}
    extractor::FeatureSubsetExtractor
    processor::DataPreprocessor
    action_clamper::DataClamper # contains action vector and clamping range

    # preallocated memory
    ϕ::Vector{Float64} # [n_predictors+1] (last one is always 1.0)

    function DynamicForestBehavior(
        forest::Ensemble,
        n_predictors::Int,
        params::DF_TrainParams,
        pre::DF_PreallocatedData,
        )

        extractor = pre.extractor
        action_clamper = pre.action_clamper
        processor = params.use_PCA ? pre.preprocess_pca : pre.preprocess

        retval = new()
        retval.forest = forest
        retval.targets = params.targets

        # copy the processor and ensure that the copied io arrays match up
        retval.extractor = FeatureSubsetExtractor(deepcopy(extractor.x), extractor.indicators)
        retval.processor = deepcopy(processor, retval.extractor)
        retval.action_clamper = DataClamper(Array(Float64, 2),
                                   deepcopy(action_clamper.f_lo),
                                   deepcopy(action_clamper.f_hi))

        retval.ϕ = Array(Float64, n_predictors+1)

        retval
    end
end
AutomotiveDrivingModels.get_targets(behavior::DynamicForestBehavior) = behavior.targets
AutomotiveDrivingModels.get_indicators(behavior::DynamicForestBehavior) = behavior.extractor.indicators

function preallocate_learning_data(
    dset::ModelTrainingData2,
    params::DF_TrainParams)

    DF_PreallocatedData(dset, params)
end

function _condition_predictor_mean!(
    normal::MvNormal,
    predictor_indeces::Vector{Int},
    A::Matrix{Float64},
    behavior::DynamicForestBehavior,
    )

    ϕ = behavior.ϕ

    for (i,p) in enumerate(predictor_indeces)
        ϕ[i] = behavior.processor.z[p]
    end
    ϕ[end] = 1.0

    normal.μ[:] = A * ϕ

    nothing
end
function select_action(
    behavior::DynamicForestBehavior,
    runlog::RunLog,
    sn::StreetNetwork,
    colset::UInt,
    frame::Int
    )

    observe!(behavior.extractor, runlog, sn, colset, frame)
    process!(behavior.processor)

    # pick a tree at random
    tree_index = rand(1:length(behavior.forest.trees))

    # sample from the MvNorm for said tree
    leaf = apply_tree(behavior.forest.trees[tree_index], behavior.processor.z)::AutoregressiveMvNormLeaf
    _condition_predictor_mean!(leaf.m, leaf.predictor_indeces, leaf.A, behavior)

    Distributions._rand!(leaf.m, behavior.action_clamper.x)
    process!(behavior.action_clamper)

    action_lat = behavior.action_clamper.x[1]
    action_lon = behavior.action_clamper.x[2]

    (action_lat, action_lon)
end

function _set_and_process_action!(behavior::DynamicForestBehavior, action_lat::Float64, action_lon::Float64)
    behavior.action_clamper.x[1] = action_lat
    behavior.action_clamper.x[2] = action_lon
    process!(behavior.action_clamper)
    behavior.action_clamper
end
function _calc_action_loglikelihood(behavior::DynamicForestBehavior)

    #=
    Action should be in behavior.action_clamper and should already have been processed
    Observations should be in behavior.processor.z and should already have been processed
    =#

    ntrees = length(behavior.forest.trees)
    log_ntrees = log(ntrees)

    total_probability_density = 0.0
    logl = 0.0
    for tree in behavior.forest.trees
        leaf = apply_tree(tree, behavior.processor.z)::AutoregressiveMvNormLeaf
        _condition_predictor_mean!(leaf.m, leaf.predictor_indeces, leaf.A, behavior)
        normal  = leaf.m

        p = pdf(normal, behavior.action_clamper.x)/ntrees

        if total_probability_density == 0.0
            logl = logpdf(normal, behavior.action_clamper.x) - log_ntrees
        elseif p < total_probability_density
            # NOTE: log(a + b) = log a + log(1 + b/a)
            logl = logl + log1p(p/total_probability_density)
        else
            # NOTE: log(a + b) = log b + log(1 + a/b)
            logl = (logpdf(normal, behavior.action_clamper.x) - log_ntrees) + log1p(total_probability_density/p)
        end
        total_probability_density += p
    end

    logl
end

function calc_action_loglikelihood(
    behavior::DynamicForestBehavior,
    runlog::RunLog,
    sn::StreetNetwork,
    colset::UInt,
    frame::Int,
    action_lat::Float64,
    action_lon::Float64,
    )

    #=
    Compute the log-likelihood of the action taken during a single frame
    given the VehicleBehaviorGaussian.
    =#

    observe!(behavior.extractor, runlog, sn, colset, frame)
    process!(behavior.processor)

    _set_and_process_action!(behavior, action_lat, action_lon)

    _calc_action_loglikelihood(behavior)
end
function calc_action_loglikelihood(
    behavior::DynamicForestBehavior,
    features::DataFrame,
    frameind::Integer,
    )

    observe!(behavior.extractor, features, frameind)
    process!(behavior.processor)

    action_lat = features[frameind, symbol(behavior.targets.lat)]::Float64
    action_lon = features[frameind, symbol(behavior.targets.lon)]::Float64
    _set_and_process_action!(behavior, action_lat, action_lon)

    _calc_action_loglikelihood(behavior)
end

function train(
    training_data::ModelTrainingData2,
    preallocated_data::DF_PreallocatedData,
    params::DF_TrainParams,
    foldset::FoldSet,
    )

    extractor = preallocated_data.extractor

    ntrees = params.ntrees
    max_tree_depth = params.max_tree_depth
    n_split_tries = params.n_split_tries
    min_samples_split = params.min_samples_split
    min_samples_leaves = params.min_samples_leaves
    n_random_predictor_samples = params.n_random_predictor_samples
    n_autoregression_predictors = params.n_autoregression_predictors
    min_split_improvement = params.min_split_improvement
    partial_sampling = params.partial_sampling
    autogression_coef = params.autogression_coef

    X_orig = params.use_PCA ? preallocated_data.Z : preallocated_data.X
    Y_orig = preallocated_data.Y

    # TODO(tim): do this without globals
    global DEFAULT_AUTOREGRESSION_CONSTANT = autogression_coef
    global DEFAULT_NUM_AUTOREGRESSION_PREDICTORS = n_autoregression_predictors
    global DEFAULT_NUM_PREDICTOR_SAMPLES = n_random_predictor_samples
    global GLOBAL_DIAG_SHRINKAGE_THRESOLD = params.diag_shrinkage_threshold

    build_tree_params = BuildTreeParameters(
        round(Int, sqrt(length(extractor.indicators))),
        max_tree_depth,
        min_samples_split,
        min_samples_leaves,
        min_split_improvement,
        n_split_tries,
        LossFunction_MSE,
        AutoregressiveMvNormLeaf
        )

    nframes = length(foldset)
    X = Array(Float64, nframes, size(X_orig, 1)) # NOTE: this is transpose of X_orig
    y = Array(Float64, 2, nframes)

    i = 0
    for frame in foldset
        i += 1

        for j in 1 : size(X_orig, 1)
            X[i,j] = X_orig[j,frame]
        end

        y[1,i] = Y_orig[1,frame]
        y[2,i] = Y_orig[2,frame]
    end
    @assert(i == nframes)

    ensemble = build_forest(y, X, ntrees, build_tree_params, partial_sampling)
    DynamicForestBehavior(ensemble, n_autoregression_predictors, params, preallocated_data)
end

end # end module


text:
import MathOptInterface
import SparseArrays

const MOI = MathOptInterface

# Supported scalar sets
const SCALAR_SETS = Union{
    MOI.GreaterThan{Float64},
    MOI.LessThan{Float64},
    MOI.EqualTo{Float64},
    MOI.Interval{Float64}
}

# Maps Clp's parameters to getter/setter function
const CLP_OPTION_MAP = Dict(
    :PrimalTolerance => (Clp_primalTolerance, Clp_setPrimalTolerance),
    :DualTolerance => (Clp_dualTolerance, Clp_setDualTolerance),
    :DualObjectiveLimit => (Clp_dualObjectiveLimit, Clp_setDualObjectiveLimit),
    :MaximumIterations => (
        # TODO(odow): Clp doesn't follow convention with maximumIterations.
        maximumIterations, Clp_setMaximumIterations
    ),
    :MaximumSeconds => (Clp_maximumSeconds, Clp_setMaximumSeconds),
    :LogLevel => (Clp_logLevel, Clp_setLogLevel),
    :Scaling => (Clp_scalingFlag, Clp_scaling),
    :Perturbation => (Clp_perturbation, Clp_setPerturbation),
    :Algorithm => (Clp_algorithm, Clp_setAlgorithm)
)

const SOLVE_OPTION_MAP = Dict(
   :PresolveType => (ClpSolve_getPresolveType, ClpSolve_setPresolveType),
   :SolveType => (ClpSolve_getSolveType, ClpSolve_setSolveType),
   :InfeasibleReturn => (ClpSolve_infeasibleReturn, ClpSolve_setInfeasibleReturn)
)

mutable struct Optimizer <: MOI.AbstractOptimizer
    inner::Ptr{Cvoid}
    solver_options::Ptr{Cvoid}
    options_set::Set{Symbol}
    optimize_called::Bool
    solve_time::Float64
    # Work-around for upstream bug in Clp:
    maximumSeconds::Float64

    function Optimizer(; kwargs...)
        model = new(
            Clp_newModel(),
            ClpSolve_new(),
            Set{Symbol}(),
            false,
            0.0,
            -1.0,
        )
        for (key, value) in kwargs
            MOI.set(model, MOI.RawParameter(String(key)), value)
        end
        finalizer(model) do m
            Clp_deleteModel(m.inner)
            ClpSolve_delete(m.solver_options)
        end
        return model
    end
end

# ====================
#   empty functions
# ====================

function MOI.is_empty(model::Optimizer)
    # A problem is empty if it has no variable and no linear constraints
    return (Clp_getNumRows(model.inner) == 0) && (Clp_getNumCols(model.inner) == 0)
end

function MOI.empty!(model::Optimizer)
    old_model = model.inner
    model.inner = Clp_newModel()
    model.optimize_called = false
    model.solve_time = 0.0
    # Copy parameters from old model into new model
    for key in model.options_set
        getter, setter = CLP_OPTION_MAP[key]
        setter(model.inner, getter(old_model))
    end
    Clp_deleteModel(old_model)
    # Work-around for maximumSeconds
    Clp_setMaximumSeconds(model.inner, model.maximumSeconds)
    return
end

MOI.get(::Optimizer, ::MOI.SolverName) = "Clp"

# TODO: improve type-stability of `MOI.RawParameter`-related methods.
MOI.supports(::Optimizer, param::MOI.RawParameter) = true

function MOI.set(model::Optimizer, param::MOI.RawParameter, value)
    key = Symbol(param.name)

    # Setters for `SolveType`, `PresolveType` and `InfeasibleReturn`
    # are handled separately, see https://github.com/jump-dev/Clp.jl/issues/91
    if haskey(CLP_OPTION_MAP, key)
        push!(model.options_set, key)
        CLP_OPTION_MAP[key][2](model.inner, value)
    elseif key == :SolveType
        ClpSolve_setSolveType(model.solver_options, value, -1)
    elseif key == :PresolveType
        ClpSolve_setPresolveType(model.solver_options, value, -1)
    elseif key == :InfeasibleReturn
        ClpSolve_setInfeasibleReturn(model.solver_options, value)
    else
        throw(MOI.UnsupportedAttribute(param))
    end
    return
end

function MOI.get(model::Optimizer, param::MOI.RawParameter)
    key = Symbol(param.name)
    if haskey(CLP_OPTION_MAP, key)
        return CLP_OPTION_MAP[key][1](model.inner)
    elseif haskey(SOLVE_OPTION_MAP, key)
        return SOLVE_OPTION_MAP[key][1](model.solver_options)
    else
        throw(MOI.UnsupportedAttribute(param))
    end
end

MOI.supports(::Optimizer, ::MOI.Silent) = true
function MOI.set(model::Optimizer, ::MOI.Silent, value::Bool)
    push!(model.options_set, :LogLevel)
    Clp_setLogLevel(model.inner, value ? 0 : 1)
    return
end

function MOI.get(model::Optimizer, ::MOI.Silent)
    return Clp_logLevel(model.inner) == 0
end

MOI.supports(::Optimizer, ::MOI.TimeLimitSec) = true

function MOI.set(model::Optimizer, ::MOI.TimeLimitSec, value)
    push!(model.options_set, :MaximumSeconds)
    value = value === nothing ? -1.0 : value
    Clp_setMaximumSeconds(model.inner, value)
    model.maximumSeconds = value
    return
end

function MOI.get(model::Optimizer, ::MOI.TimeLimitSec)
    # TODO(odow): replace with `Clp_maximumSeconds(model.inner)` when upstream
    # is fixed.
    return model.maximumSeconds
end

MOI.supports(::Optimizer, ::MOI.NumberOfThreads) = false

# ========================================
#   Supported constraints and objectives
# ========================================

function MOI.supports_constraint(
    ::Optimizer, ::Type{MOI.ScalarAffineFunction{Float64}}, ::Type{<:SCALAR_SETS}
)
    return true
end

function MOI.supports_constraint(
    ::Optimizer, ::Type{MOI.SingleVariable}, ::Type{<:SCALAR_SETS}
)
    return true
end

MOI.supports(::Optimizer, ::MOI.ObjectiveSense) = true

function MOI.supports(
    ::Optimizer, ::MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}
)
    return true
end

# =======================
#   `copy_to` function
# =======================

_add_bounds(::Vector{Float64}, ub, i, s::MOI.LessThan{Float64}) = ub[i] = s.upper
_add_bounds(lb, ::Vector{Float64}, i, s::MOI.GreaterThan{Float64}) = lb[i] = s.lower
_add_bounds(lb, ub, i, s::MOI.EqualTo{Float64}) = lb[i], ub[i] = s.value, s.value
_add_bounds(lb, ub, i, s::MOI.Interval{Float64}) = lb[i], ub[i] = s.lower, s.upper

function _extract_bound_data(src, mapping, lb, ub, ::Type{S}) where S
    for con_index in MOI.get(
        src, MOI.ListOfConstraintIndices{MOI.SingleVariable, S}()
    )
        f = MOI.get(src, MOI.ConstraintFunction(), con_index)
        s = MOI.get(src, MOI.ConstraintSet(), con_index)
        column = mapping.varmap[f.variable].value
        _add_bounds(lb, ub, column, s)
        mapping.conmap[con_index] = MOI.ConstraintIndex{MOI.SingleVariable, S}(column)
    end
end

function _copy_to_columns(dest, src, mapping)
    x_src = MOI.get(src, MOI.ListOfVariableIndices())
    N = Cint(length(x_src))
    for i = 1:N
        mapping.varmap[x_src[i]] = MOI.VariableIndex(i)
    end

    fobj = MOI.get(src, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}())
    c = fill(0.0, N)
    for term in fobj.terms
        i = mapping.varmap[term.variable_index].value
        c[i] += term.coefficient
    end
    # Clp seems to negates the objective offset
    Clp_setObjectiveOffset(dest.inner, -fobj.constant)
    return N, c
end

_bounds(s::MOI.GreaterThan{Float64}) = (s.lower, Inf)
_bounds(s::MOI.LessThan{Float64}) = (-Inf, s.upper)
_bounds(s::MOI.EqualTo{Float64}) = (s.value, s.value)
_bounds(s::MOI.Interval{Float64}) = (s.lower, s.upper)

function add_sizehint!(vec, n)
    len = length(vec)
    return sizehint!(vec, len + n)
end

function _extract_row_data(src, mapping, lb, ub, I, J, V, ::Type{S}) where S
    row = length(I) == 0 ? 1 : I[end] + 1
    list = MOI.get(
        src, MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{Float64}, S}()
    )
    add_sizehint!(lb, length(list))
    add_sizehint!(ub, length(list))
    n_terms = 0
    fs = Array{MOI.ScalarAffineFunction{Float64}}(undef, length(list))
    for (i,c_index) in enumerate(list)
        f = MOI.get(src, MOI.ConstraintFunction(), c_index)
        fs[i] = f
        l, u = _bounds(MOI.get(src, MOI.ConstraintSet(), c_index))
        push!(lb, l - f.constant)
        push!(ub, u - f.constant)
        n_terms += length(f.terms)
    end
    add_sizehint!(I, n_terms)
    add_sizehint!(J, n_terms)
    add_sizehint!(V, n_terms)
    for (i,c_index) in enumerate(list)
        f = fs[i]#MOI.get(src, MOI.ConstraintFunction(), c_index)
        for term in f.terms
            push!(I, row)
            push!(J, Cint(mapping.varmap[term.variable_index].value))
            push!(V, term.coefficient)
        end
        mapping.conmap[c_index] = MOI.ConstraintIndex{
            MOI.ScalarAffineFunction{Float64}, S
        }(length(ub))
        row += 1
    end
    return
end

function test_data(src, dest)
    for (F, S) in MOI.get(src, MOI.ListOfConstraints())
        if !MOI.supports_constraint(dest, F, S)
            throw(MOI.UnsupportedConstraint{F, S}("Clp.Optimizer does not support constraints of type $F-in-$S."))
        end
    end
    fobj_type = MOI.get(src, MOI.ObjectiveFunctionType())
    if !MOI.supports(dest, MOI.ObjectiveFunction{fobj_type}())
        throw(MOI.UnsupportedAttribute(MOI.ObjectiveFunction(fobj_type)))
    end
end

function MOI.copy_to(
    dest::Optimizer,
    src::MOI.ModelLike;
    copy_names::Bool = false
)
    @assert MOI.is_empty(dest)
    test_data(src, dest)

    mapping = MOI.Utilities.IndexMap()
    N, c = _copy_to_columns(dest, src, mapping)
    cl, cu = fill(-Inf, N), fill(Inf, N)
    rl, ru, I, J, V = Float64[], Float64[], Cint[], Cint[], Float64[]

    _extract_bound_data(src, mapping, cl, cu, MOI.GreaterThan{Float64})
    _extract_row_data(src, mapping, rl, ru, I, J, V, MOI.GreaterThan{Float64})
    _extract_bound_data(src, mapping, cl, cu, MOI.LessThan{Float64})
    _extract_row_data(src, mapping, rl, ru, I, J, V, MOI.LessThan{Float64})
    _extract_bound_data(src, mapping, cl, cu, MOI.EqualTo{Float64})
    _extract_row_data(src, mapping, rl, ru, I, J, V, MOI.EqualTo{Float64})
    _extract_bound_data(src, mapping, cl, cu, MOI.Interval{Float64})
    _extract_row_data(src, mapping, rl, ru, I, J, V, MOI.Interval{Float64})

    M = Cint(length(rl))
    A = SparseArrays.sparse(I, J, V, M, N)
    Clp_loadProblem(
        dest.inner,
        A.n,
        A.m,
        A.colptr .- Cint(1),
        A.rowval .- Cint(1),
        A.nzval,
        cl,
        cu,
        c,
        rl,
        ru
    )

    sense = MOI.get(src, MOI.ObjectiveSense())
    if sense == MOI.MIN_SENSE
        Clp_setObjSense(dest.inner, 1)
    elseif sense == MOI.MAX_SENSE
        Clp_setObjSense(dest.inner, -1)
    else
        @assert sense == MOI.FEASIBILITY_SENSE
        Clp_setObjSense(dest.inner, 0)
    end
    return mapping
end

# ===============================
#   Optimize and post-optimize
# ===============================

function MOI.optimize!(model::Optimizer)
    t = time()
    Clp_initialSolveWithOptions(model.inner, model.solver_options)
    model.solve_time = time() - t
    model.optimize_called = true
    return
end

function MOI.get(model::Optimizer, ::MOI.SolveTime)
    return model.solve_time
end

function MOI.get(model::Optimizer, ::MOI.NumberOfVariables)
    return Clp_getNumCols(model.inner)
end

function MOI.get(model::Optimizer, attr::MOI.ObjectiveValue)
    MOI.check_result_index_bounds(model, attr)
    return Clp_getObjValue(model.inner)
end

function MOI.get(model::Optimizer, ::MOI.TerminationStatus)
    if !model.optimize_called
        return MOI.OPTIMIZE_NOT_CALLED
    end
    st = Clp_status(model.inner)
    if st == 0
        return MOI.OPTIMAL
    elseif st == 1
        return MOI.INFEASIBLE
    elseif st == 2
        return MOI.DUAL_INFEASIBLE
    elseif st == 3
        # No more granular information that "some limit is reached"
        return MOI.OTHER_LIMIT
    else
        return MOI.OTHER_ERROR
    end
end

function MOI.get(model::Optimizer, ::MOI.RawStatusString)
    if !model.optimize_called
        return "MOI.OPTIMIZE_NOT_CALLED"
    end
    st = Clp_status(model.inner)
    if st == 0
        return "0 - optimal"
    elseif st == 1
        return "1 - primal infeasible"
    elseif st == 2
        return "2 - dual infeasible"
    elseif st == 3
        return "3 - stopped on iterations etc"
    elseif st == 4
        return "4 - stopped due to errors"
    else
        error("Expected integer in [0, 4] but got $st")
    end
end

function MOI.get(model::Optimizer, ::MOI.ResultCount)
    if Clp_primalFeasible(model.inner) != 0
        return 1
    elseif Clp_dualFeasible(model.inner) != 0
        return 1
    elseif Clp_isProvenPrimalInfeasible(model.inner) != 0
        return 1
    elseif Clp_isProvenDualInfeasible(model.inner) != 0
        return 1
    end
    return 0
end

function MOI.get(model::Optimizer, attr::MOI.PrimalStatus)
    if attr.N != 1
        return MOI.NO_SOLUTION
    elseif Clp_isProvenDualInfeasible(model.inner) != 0
        return MOI.INFEASIBILITY_CERTIFICATE
    elseif Clp_primalFeasible(model.inner) != 0
        return MOI.FEASIBLE_POINT
    else
        return MOI.UNKNOWN_RESULT_STATUS
    end
end

function MOI.get(model::Optimizer, attr::MOI.DualStatus)
    if attr.N != 1
        return MOI.NO_SOLUTION
    elseif Clp_isProvenPrimalInfeasible(model.inner) != 0
        return MOI.INFEASIBILITY_CERTIFICATE
    elseif Clp_dualFeasible(model.inner) != 0
        return MOI.FEASIBLE_POINT
    else
        return MOI.UNKNOWN_RESULT_STATUS
    end
end

# ===================
#   Primal solution
# ===================

function _unsafe_wrap_clp_array(
    model::Ptr{Cvoid},
    f::Function,
    n::Integer,
    indices;
    own::Bool = false
)
    p = f(model)
    if p == C_NULL
        return map(x -> NaN, indices)
    end
    x = unsafe_wrap(Array, p, (n,); own = own)
    return x[indices]
end

function MOI.get(
    model::Optimizer, attr::MOI.VariablePrimal, x::MOI.VariableIndex
)
    MOI.check_result_index_bounds(model, attr)
    primal_status = MOI.get(model, MOI.PrimalStatus())
    if primal_status == MOI.INFEASIBILITY_CERTIFICATE
        # We claim ownership of the pointer returned by Clp_unboundedRay.
        return _unsafe_wrap_clp_array(
            model.inner,
            Clp_unboundedRay,
            Clp_getNumCols(model.inner),
            x.value;
            own = true,
        )
    elseif primal_status == MOI.FEASIBLE_POINT
        return _unsafe_wrap_clp_array(
            model.inner,
            Clp_getColSolution,
            Clp_getNumCols(model.inner),
            x.value,
        )
    else
        error("Primal solution not available")
    end
end

function MOI.get(
    model::Optimizer, attr::MOI.VariablePrimal, xs::Vector{MOI.VariableIndex}
)
    MOI.check_result_index_bounds(model, attr)
    col_indices = [idx.value for idx in xs]
    primal_status = MOI.get(model, MOI.PrimalStatus())
    if primal_status == MOI.INFEASIBILITY_CERTIFICATE
        # We claim ownership of the pointer returned by Clp_unboundedRay.
        return _unsafe_wrap_clp_array(
            model.inner,
            Clp_unboundedRay,
            Clp_getNumCols(model.inner),
            col_indices;
            own = true,
        )
    elseif primal_status == MOI.FEASIBLE_POINT
        return _unsafe_wrap_clp_array(
            model.inner,
            Clp_getColSolution,
            Clp_getNumCols(model.inner),
            col_indices,
        )
    else
        error("Primal solution not available")
    end
end

# TODO: What happens if model is unbounded / infeasible?
function MOI.get(
    model::Optimizer,
    attr::MOI.ConstraintPrimal,
    c::MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64}, <:SCALAR_SETS}
)
    MOI.check_result_index_bounds(model, attr)
    return _unsafe_wrap_clp_array(
        model.inner,
        Clp_getRowActivity,
        Clp_getNumRows(model.inner),
        c.value,
    )
end

# TODO: What happens if model is unbounded / infeasible?
function MOI.get(
    model::Optimizer,
    attr::MOI.ConstraintPrimal,
    c::MOI.ConstraintIndex{MOI.SingleVariable, <:SCALAR_SETS}
)
    MOI.check_result_index_bounds(model, attr)
    return MOI.get(model, MOI.VariablePrimal(), MOI.VariableIndex(c.value))
end

# =================
#   Dual solution
# =================

# If sense is maximize, we negate all the duals to follow MOI conventions
# Feasibility problems are treated as a minimization

function MOI.get(
    model::Optimizer,
    attr::MOI.ConstraintDual,
    c::MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64}, <:SCALAR_SETS}
)
    MOI.check_result_index_bounds(model, attr)
    sense = (Clp_getObjSense(model.inner) == -1) ? -1 : 1
    dual_status = MOI.get(model, MOI.DualStatus())
    if dual_status == MOI.FEASIBLE_POINT
        dsol = _unsafe_wrap_clp_array(
            model.inner,
            Clp_getRowPrice,
            Clp_getNumRows(model.inner),
            c.value,
        )
        return sense * dsol
    elseif dual_status == MOI.INFEASIBILITY_CERTIFICATE
        # We claim ownership of the pointer returned by Clp_infeasibilityRay.
        dsol = _unsafe_wrap_clp_array(
            model.inner,
            Clp_infeasibilityRay,
            Clp_getNumRows(model.inner),
            c.value;
            own = true,
        )
        return -sense * dsol
    else
        error("Dual solution not available")
    end
end

# TODO: what happens if problem is unbounded / infeasible?
function MOI.get(
    model::Optimizer,
    attr::MOI.ConstraintDual,
    c::MOI.ConstraintIndex{MOI.SingleVariable, MOI.LessThan{Float64}}
)
    MOI.check_result_index_bounds(model, attr)
    rc = _unsafe_wrap_clp_array(
        model.inner,
        Clp_getReducedCost,
        Clp_getNumCols(model.inner),
        c.value,
    )
    sense = (Clp_getObjSense(model.inner) == -1) ? -1 : 1
    if sense == 1 && rc <= 0.0
        return rc
    elseif sense == -1 && rc >= 0.0
        return -rc
    else
        return 0.0
    end
end

function MOI.get(
    model::Optimizer,
    attr::MOI.ConstraintDual,
    c::MOI.ConstraintIndex{MOI.SingleVariable, MOI.GreaterThan{Float64}}
)
    MOI.check_result_index_bounds(model, attr)
    rc = _unsafe_wrap_clp_array(
        model.inner,
        Clp_getReducedCost,
        Clp_getNumCols(model.inner),
        c.value,
    )
    sense = (Clp_getObjSense(model.inner) == -1) ? -1 : 1
    if sense == 1 && rc >= 0.0
        return rc
    elseif sense == -1 && rc <= 0.0
        return -rc
    else
        return 0.0
    end
end

function MOI.get(
    model::Optimizer,
    attr::MOI.ConstraintDual,
    c::MOI.ConstraintIndex{
        MOI.SingleVariable,
        <:Union{MOI.Interval{Float64}, MOI.EqualTo{Float64}}
    }
)
    MOI.check_result_index_bounds(model, attr)
    sense = (Clp_getObjSense(model.inner) == -1) ? -1 : 1
    rc = _unsafe_wrap_clp_array(
        model.inner,
        Clp_getReducedCost,
        Clp_getNumCols(model.inner),
        c.value,
    )
    return sense * rc
end


text:
function moving_object(object; start_point, final_point, frames = 1:100)
    act!(object, Action(frames, polyout(2), anim_translate(start_point, final_point)))
end

function moving_text(;
    start_point,
    final_point,
    frames = 1:100,
    body = "stdtxt",
    size = 200,
    valign = :baseline,
    halign = :center,
    kwargs...,
)

    ob = Object(
        frames,
        (args...) -> draw_text(
            args...,
            body = body,
            size = size,
            pos = start_point,
            valign = valign,
            halign = halign,
        ),
    )

    act!(ob, Action(1:length(frames), polyout(2), anim_translate(start_point, final_point)))
    ob
end


text:
### Abstract SMMALA state

abstract SMMALAState{F<:VariateForm} <: LMCSamplerState{F}

### SMMALA state subtypes

## UnvSMMALAState holds the internal state ("local variables") of the SMMALA sampler for univariate parameters

type UnvSMMALAState <: SMMALAState{Univariate}
  pstate::ParameterState{Continuous, Univariate} # Parameter state used internally by SMMALA
  tune::MCTunerState
  sqrttunestep::Real
  ratio::Real
  μ::Real
  newinvtensor::Real
  oldinvtensor::Real
  cholinvtensor::Real
  newfirstterm::Real
  oldfirstterm::Real

  function UnvSMMALAState(
    pstate::ParameterState{Continuous, Univariate},
    tune::MCTunerState,
    sqrttunestep::Real,
    ratio::Real,
    μ::Real,
    newinvtensor::Real,
    oldinvtensor::Real,
    cholinvtensor::Real,
    newfirstterm::Real,
    oldfirstterm::Real
  )
    if !isnan(ratio)
      @assert ratio > 0 "Acceptance ratio should be positive"
    end
    if !isnan(sqrttunestep)
      @assert sqrttunestep > 0 "Square root of tuned drift step is not positive"
    end
    new(pstate, tune, sqrttunestep, ratio, μ, newinvtensor, oldinvtensor, cholinvtensor, newfirstterm, oldfirstterm)
  end
end

UnvSMMALAState(pstate::ParameterState{Continuous, Univariate}, tune::MCTunerState=BasicMCTune()) =
  UnvSMMALAState(pstate, tune, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN)

## MuvSMMALAState holds the internal state ("local variables") of the SMMALA sampler for multivariate parameters

type MuvSMMALAState <: SMMALAState{Multivariate}
  pstate::ParameterState{Continuous, Multivariate} # Parameter state used internally by SMMALA
  tune::MCTunerState
  sqrttunestep::Real
  ratio::Real
  μ::RealVector
  newinvtensor::RealMatrix
  oldinvtensor::RealMatrix
  cholinvtensor::RealLowerTriangular
  newfirstterm::RealVector
  oldfirstterm::RealVector

  function MuvSMMALAState(
    pstate::ParameterState{Continuous, Multivariate},
    tune::MCTunerState,
    sqrttunestep::Real,
    ratio::Real,
    μ::RealVector,
    newinvtensor::RealMatrix,
    oldinvtensor::RealMatrix,
    cholinvtensor::RealLowerTriangular,
    newfirstterm::RealVector,
    oldfirstterm::RealVector
  )
    if !isnan(ratio)
      @assert ratio > 0 "Acceptance ratio should be positive"
    end
    if !isnan(sqrttunestep)
      @assert sqrttunestep > 0 "Square root of tuned drift step is not positive"
    end
    new(pstate, tune, sqrttunestep, ratio, μ, newinvtensor, oldinvtensor, cholinvtensor, newfirstterm, oldfirstterm)
  end
end

MuvSMMALAState(pstate::ParameterState{Continuous, Multivariate}, tune::MCTunerState=BasicMCTune()) =
  MuvSMMALAState(
    pstate,
    tune,
    NaN,
    NaN,
    Array(eltype(pstate), pstate.size),
    Array(eltype(pstate), pstate.size, pstate.size),
    Array(eltype(pstate), pstate.size, pstate.size),
    RealLowerTriangular(Array(eltype(pstate), pstate.size, pstate.size)),
    Array(eltype(pstate), pstate.size),
    Array(eltype(pstate), pstate.size)
  )

### Metropolis-adjusted Langevin Algorithm (SMMALA)

immutable SMMALA <: LMCSampler
  driftstep::Real
  transform::Union{Function, Void} # Function for transforming metric tensor to a positive-definite matrix

  function SMMALA(driftstep::Real, transform::Union{Function, Void})
    @assert driftstep > 0 "Drift step is not positive"
    new(driftstep, transform)
  end
end

SMMALA(driftstep::Real=1.) = SMMALA(driftstep, nothing)

### Initialize SMMALA sampler

## Initialize parameter state

function initialize!(
  pstate::ParameterState{Continuous, Univariate},
  parameter::Parameter{Continuous, Univariate},
  sampler::SMMALA,
  outopts::Dict
)
  parameter.uptotensorlogtarget!(pstate)
  @assert isfinite(pstate.logtarget) "Log-target not finite: initial value out of support"
  @assert all(isfinite(pstate.gradlogtarget)) "Gradient of log-target not finite: initial values out of support"
  @assert all(isfinite(pstate.tensorlogtarget)) "Tensor of log-target not finite: initial values out of support"

  if !isempty(outopts[:diagnostics])
    pstate.diagnostickeys = copy(outopts[:diagnostics])
    pstate.diagnosticvalues = Array(Any, length(pstate.diagnostickeys))
  end
end

function initialize!(
  pstate::ParameterState{Continuous, Multivariate},
  parameter::Parameter{Continuous, Multivariate},
  sampler::SMMALA,
  outopts::Dict
)
  parameter.uptotensorlogtarget!(pstate)
  if sampler.transform != nothing
    pstate.tensorlogtarget[:, :] = sampler.transform(pstate.tensorlogtarget)
  end
  @assert isfinite(pstate.logtarget) "Log-target not finite: initial value out of support"
  @assert all(isfinite(pstate.gradlogtarget)) "Gradient of log-target not finite: initial values out of support"
  @assert all(isfinite(pstate.tensorlogtarget)) "Tensor of log-target not finite: initial values out of support"

  if !isempty(outopts[:diagnostics])
    pstate.diagnostickeys = copy(outopts[:diagnostics])
    pstate.diagnosticvalues = Array(Any, length(pstate.diagnostickeys))
  end
end

## Initialize SMMALA state

function sampler_state(
  parameter::Parameter{Continuous, Univariate},
  sampler::SMMALA,
  tuner::MCTuner,
  pstate::ParameterState{Continuous, Univariate},
  vstate::VariableStateVector
)
  sstate = UnvSMMALAState(generate_empty(pstate), tuner_state(parameter, sampler, tuner))
  sstate.sqrttunestep = sqrt(sstate.tune.step)
  sstate.oldinvtensor = inv(pstate.tensorlogtarget)
  sstate.cholinvtensor = chol(sstate.oldinvtensor, :L)
  sstate.oldfirstterm = sstate.oldinvtensor*pstate.gradlogtarget
  sstate
end

function sampler_state(
  parameter::Parameter{Continuous, Multivariate},
  sampler::SMMALA,
  tuner::MCTuner,
  pstate::ParameterState{Continuous, Multivariate},
  vstate::VariableStateVector
)
  sstate = MuvSMMALAState(generate_empty(pstate), tuner_state(parameter, sampler, tuner))
  sstate.sqrttunestep = sqrt(sstate.tune.step)
  sstate.oldinvtensor[:, :] = inv(pstate.tensorlogtarget)
  sstate.cholinvtensor[:, :] = chol(sstate.oldinvtensor, Val{:L})
  sstate.oldfirstterm[:] = sstate.oldinvtensor*pstate.gradlogtarget
  sstate
end

### Reset SMMALA sampler

## Reset parameter state

function reset!(
  pstate::ParameterState{Continuous, Univariate},
  x::Real,
  parameter::Parameter{Continuous, Univariate},
  sampler::SMMALA
)
  pstate.value = x
  parameter.uptotensorlogtarget!(pstate)
end

function reset!(
  pstate::ParameterState{Continuous, Multivariate},
  x::RealVector,
  parameter::Parameter{Continuous, Multivariate},
  sampler::SMMALA
)
  pstate.value = copy(x)
  parameter.uptotensorlogtarget!(pstate)
end

## Reset sampler state

function reset!(
  sstate::UnvSMMALAState,
  pstate::ParameterState{Continuous, Univariate},
  parameter::Parameter{Continuous, Univariate},
  sampler::SMMALA,
  tuner::MCTuner
)
  reset!(sstate.tune, sampler, tuner)
  sstate.sqrttunestep = sqrt(sstate.tune.step)
  sstate.oldinvtensor = inv(pstate.tensorlogtarget)
  sstate.cholinvtensor = chol(sstate.oldinvtensor, :L)
  sstate.oldfirstterm = sstate.oldinvtensor*pstate.gradlogtarget
end

function reset!(
  sstate::MuvSMMALAState,
  pstate::ParameterState{Continuous, Multivariate},
  parameter::Parameter{Continuous, Multivariate},
  sampler::SMMALA,
  tuner::MCTuner
)
  reset!(sstate.tune, sampler, tuner)
  sstate.sqrttunestep = sqrt(sstate.tune.step)
  sstate.oldinvtensor[:, :] = inv(pstate.tensorlogtarget)
  sstate.cholinvtensor[:, :] = chol(sstate.oldinvtensor, Val{:L})
  sstate.oldfirstterm[:] = sstate.oldinvtensor*pstate.gradlogtarget
end

Base.show(io::IO, sampler::SMMALA) = print(io, "SMMALA sampler: drift step = $(sampler.driftstep)")

Base.writemime(io::IO, ::MIME"text/plain", sampler::SMMALA) = show(io, sampler)


text:
using ClimateMachine.ODESolvers
using ClimateMachine.SystemSolvers

const slow_mrrk_methods =
    ((LSRK54CarpenterKennedy, 4), (LSRK144NiegemannDiehlBusch, 4))
const fast_mrrk_methods = (
    (LSRK54CarpenterKennedy, 4),
    (LSRK144NiegemannDiehlBusch, 4),
    (SSPRK33ShuOsher, 3),
    (SSPRK34SpiteriRuuth, 3),
)
const explicit_methods = (
    (LSRK54CarpenterKennedy, 4),
    (LSRK144NiegemannDiehlBusch, 4),
    (LS3NRK44Classic, 4),
    (LS3NRK33Heuns, 3),
    (SSPRK22Heuns, 2),
    (SSPRK22Ralstons, 2),
    (SSPRK33ShuOsher, 3),
    (SSPRK34SpiteriRuuth, 3),
    (LSRKEulerMethod, 1),
)

const imex_methods_lowstorage_compatible = (
    # Low-storage variant methods have an assumption that the
    # explicit and implicit rhs/time-scaling coefficients (B/C vectors)
    # in the Butcher tables are the same.
    (ARK1ForwardBackwardEuler, 1),
    (ARK2ImplicitExplicitMidpoint, 2),
    (ARK2GiraldoKellyConstantinescu, 2),
    (ARK437L2SA1KennedyCarpenter, 4),
    (ARK548L2SA2KennedyCarpenter, 5),
    (DBM453VoglEtAl, 3),
)
const imex_methods_naivestorage_compatible = (
    imex_methods_lowstorage_compatible...,
    # Some methods can only be used with the `NaiveVariant` storage
    # scheme since, in general, ARK methods can have different time-scaling/rhs-scaling
    # coefficients (C/B vectors in the Butcher tables). For future reference,
    # any other ARK-type methods that have more general Butcher tables
    # (but with same number of stages) should be tested here:
    (Trap2LockWoodWeller, 2),
)

const mis_methods =
    ((MIS2, 2), (MIS3C, 2), (MIS4, 3), (MIS4a, 3), (TVDMISA, 2), (TVDMISB, 2))

const mrigark_erk_methods = ((MRIGARKERK33aSandu, 3), (MRIGARKERK45aSandu, 4))

const mrigark_irk_methods = (
    (MRIGARKESDIRK24LSA, 2),
    (MRIGARKESDIRK23LSA, 2),
    (MRIGARKIRK21aSandu, 2),
    (MRIGARKESDIRK34aSandu, 3),
    (MRIGARKESDIRK46aSandu, 4),
)

const fast_mrigark_methods =
    ((LSRK54CarpenterKennedy, 4), (LSRK144NiegemannDiehlBusch, 4))

struct DivideLinearSolver <: AbstractSystemSolver end
function SystemSolvers.prefactorize(
    linearoperator!,
    ::DivideLinearSolver,
    args...,
)
    linearoperator!
end
function SystemSolvers.linearsolve!(
    linearoperator!,
    preconditioner,
    ::DivideLinearSolver,
    Qtt,
    Qhat,
    args...,
)
    @. Qhat = 1 / Qhat
    linearoperator!(Qtt, Qhat, args...)
    @. Qtt = 1 / Qtt
end