text:
module BeastNames

const TAXA = "taxa"
const TAXON = "taxon"
const ID = "id"
const NAME = "name"
const MISSING_VAL = "NA"
const NEWICK = "newick"
const MCMC = "mcmc"
const LOG = "log"
const FILENAME = "fileName"
const CHAINLENGTH = "chainLength"
const LOGEVERY = "logEvery"
const DATE = "date"
const OPERATORS = "operators"
const ATTR = "attr"
const MULTIVARIATE_DIFFUSION_MODEL = "multivariateDiffusionModel"
const PRECISION_MATRIX = "precisionMatrix"
const MATRIX_PARAMETER = "matrixParameter"
const PARAMETER = "parameter"
const VALUE = "value"
const SAMPLING_PRECISION = "samplingPrecision"
const REPEATED_MEASURES = "repeatedMeasuresModel"
const WEIGHT = "weight"
const BEAST = "beast"

const TRAIT_VALIDATION = "traitValidation"
const LOG_SUM = "logSum"
const TRAIT_VALIDATION_PROVIDER = "traitValidationProvider"
const CROSS_VALIDATION = "crossValidation"

const TRAIT_NAME = "traitName"
const INFERRED_TRAIT = "inferredTrait"
const USE_TREE_TRAITS = "useTreeTraits"
const TRAIT_DATA_LIKELIHOOD = "traitDataLikelihood"
const IDREF = "idref"
const TRAIT_PARAMETER = "traitParameter"
const FILELOG = "fileLog"
const MULT_DIST_LIKELIHOOD = "multivariateDistributionLikelihood"
const DISTRIBUTION = "distribution"
const MVM_DISTRIBUTION = "multivariateNormalDistributionModel"
const MEAN_PARAMETER = "meanParameter"
const PRECISION_PARAMETER = "precisionParameter"
const DATA = "data"
const LEAF_TRAIT_PARAMETER = "leafTraitParameter"
const TREE_MODEL = "treeModel"
const LIKELIHOOD = "likelihood"
const POSTERIOR = "joint"
const TRUE = "true"
const FALSE = "false"
const FAST_MATRIX_PARAMETER = "fastMatrixParameter"
const ROWS = "rows"
const COLUMNS = "columns"

const DEFAULT_TRAIT_NAME = "traits"
const DEFAULT_TREE_NAME = "startingTree"

const FIXED_TREE = "fixedTree"
const TREE = "tree"
const ROOT_HEIGHT = "rootHeight"
const INTERNAL_NODES = "internalNodes"
const ROOT_NODE = "rootNode"
const NODE_HEIGHTS = "nodeHeights"
const INTERNAL_NODE_HEIGHTS = "internalNodeHeights"
const ALL_INTENAL_NODE_HEIGHTS = "allInternalNodeHeights"
const NODE_TRAITS = "nodeTraits"
const LEAF_NODES = "leafNodes"
const TRAIT_DIM = "traitDimension"
const AS_MATRIX = "asMatrix"
const LEAF_TRAITS = "leafTraits"

const MULTIVARIATE_WISHART_PRIOR = "multivariateWishartPrior"
const SCALE_MATRIX = "scaleMatrix"
const DF = "df"
const DEFAULT_MBD_PRIOR = "precisionPrior"

const DIFFUSION_ID = "diffusionModel"
const DIAGONAL_MATRIX = "DiagonalMatrix"
const DIFF_PREC_ID = "diffusion.precision"
const DIMENSION = "dimension"
const LOWER = "lower"
const UPPER = "upper"

const DEFAULT_RM_NAME = "repeatedMeasures"
const DEFAULT_RM_PREC_NAME = "residualPrecision"
const DEFAULT_RM_PREC_PRIOR_NAME = "residualPrior"

const CACHE_BRANCHES = "cacheBranches"
const ALLOW_IDENTICAL = "allowIdentical"
const USE_TREE_LENGTH = "useTreeLength"
const SCALE_BY_TIME = "scaleByTime"
const REPORT_AS_MULTIVARIATE = "reportAsMultivariate"
const INTEGRATE_INTERNAL_TRAITS = "integrateInternalTraits"
const FACTOR_TRAIT_NAME = "factors"
const CONJUGATE_ROOT_PRIOR = "conjugateRootPrior"
const PSS = "priorSampleSize"
const TRAIT_DATA_LIKELIHOOD_ID = "traitLikelihood"
const PSS_VAL = 0.001
const ALLOW_SINGULAR = "allowSingular"
const STANDARDIZE = "standardize"

const WISHART_STATISTICS = "wishartStatistics"
const REPEATED_MEASURES_WISHART_STATISTICS = "RepeatedMeasuresWishartStatistics"
const PRECISION_GIBBS_OPERATOR = "precisionGibbsOperator"
const COMPOUND_PRECISION_OPERATOR = "compoundPrecisionOperator"
const DIFFUSION_OPERATOR = "diffusionOperator"
const RESIDUAL_OPERATOR = "residualOperator"

const PRIOR = "prior"
const PRIOR_ID = "prior"
const LIKELIHOOD_ID = "likelihood"
const SCREEN_LOGEVERY = 1000
const FILE_LOGEVERY = 1000
const SCREEN_LOG_ID = "screenLog"
const FILE_LOG_ID = "fileLog"
const OVERWRITE = "overwrite"
const AUTO_OPTIMIZE = "autoOptimize"
const DEFAULT_FILENAME = "fileLog"

const COLUMN = "column"
const LABEL = "label"
const DP = "dp"
const WIDTH = "width"
const DP_DEFAULT = 4
const WIDTH_DEFAULT = 12

const TIMER = "timer"
const REPORT = "report"
const PROPERTY = "property"

const MATRIX_INVERSE = "matrixInverse"

const DIRECTION = "direction"
const FORWARDS = "forwards"

const USING_HEIGHTS = "usingHeights"
const USING_DATES = "usingDates"
const FIX_HEIGHTS = "fixHeights"

const CORRELATION_MATRIX = "correlationMatrix"
const INVERT = "invert"

const FIRE_PARAMETER_CHANGED = "fireParameterChanged"

const VARIANCE_PROPORTION_STATISTIC = "varianceProportionStatistic"
const MATRIX_RATIO = "matrixRatio"
const COHERITABILITY = "coheritability"
const USE_EMPIRICAL_VARIANCE = "useEmpiricalVariance"

const FULL_EVALUATION = "fullEvaluation"

const MATRIX_VALIDATION = "matrixValidation"
const TRUE_PARAMETER = "trueParameter"
const INFERRED_PARAMETER = "inferredParameter"

const MODEL_EXTENSION_TRAIT_LOGGER = "modelExtensionTraitLogger"
const DIMENSIONS = "dimensions"
const ALL = "all"
const MISSING = "missing"

const INTEGRATED_FACTORS = "integratedFactorModel"
const DEFAULT_IF_NAME = "factorModel"
const LOADINGS = "loadings"
const DEFAULT_LOADINGS_ID = "L"
const PRECISION = "precision"
const FACTOR_PRECISION = "factorPrecision"

const DISTRIBUTION_LIKELIHOOD = "distributionLikelihood"
const NORMAL_DISTRIBUTION_MODEL = "normalDistributionModel"
const MEAN = "mean"
const STDEV = "stdev"

const HMC = "hamiltonianMonteCarloOperator"
const GRAD_CHECK_COUNT = "gradientCheckCount"
const JOINT_GRADIENT = "jointGradient"
const GRADIENT = "gradient"
const INTEGRATED_FACTOR_GRADIENT = "integratedFactorAnalysisLoadingsGradient"

const NSTEPS = "nSteps"
const STEP_SIZE = "stepSize"
const DRAW_VARIANCE = "drawVariance"
const GRAD_TOLERANCE = "gradientCheckTolerance"

const GAMMA_PRIOR = "gammaPrior"
const SCALE = "scale"
const SHAPE = "shape"

const NGP_OPERATOR = "normalGammaPrecisionGibbsOperator"
const NORMAL_EXTENSION = "normalExtension"
const TREE_TRAIT_NAME = "treeTraitName"

const DEFAULT_FACTOR_NAME = "factors"

const LOADINGS_GIBBS_OP = "loadingsGibbsOperator"
const RANDOM_SCAN = "randomScan"
const NEW_MODE = "newMode"
const CONSTRAINT = "constraint"
const SPARSITY = "sparsity"
const NONE = "none"
const NORMAL_PRIOR = "normalPrior"

const LATENT_FACTOR_MODEL = "latentFactorModel"
const FACTOR_NUMBER = "factorNumber"
const SCALE_DATA = "scaleData"
const FACTORS = "factors"
const DATA_FROM_TIPS = "dataAndMissingFromTreeTips"
const ROW_PRECISION = "rowPrecision"
const COL_PRECISION = "columnPrecision"
const LAT_FAC_PREC_OP = "latentFactorModelPrecisionOperator"
const FACTOR_TREE_GIBBS_OP = "factorTreeGibbsOperator"

const MULTIVARIATE_TRAIT_LIKELIHOOD = "multivariateTraitLikelihood"

const TRAIT_LOGGER = "traitLogger"
const TAXON_EXPLICIT = "taxonNameExplicit"
const NODES = "nodes"
const EXTERNAL = "external"

const PRODUCT_PARAMETER = "productParameter"
const TRANSFORMED_PARAMETER = "transformedParameter"
const TYPE = "type"
const POWER = "power"
const POWER_TRANSFORM = "powerTransform"
const TRANSFORM = "transform"

const BAYESIAN_BRIDGE = "bayesianBridge"
const GLOBAL_SCALE = "globalScale"
const EXPONENT = "exponent"
const LOCAL_SCALE = "localScale"

const MATRIX_SHRINKAGE_LIKELIHOOD = "matrixShrinkageLikelihood"
const ROW_PRIORS = "rowPriors"

const SCALE_OPERATOR = "scaleOperator"
const SCALE_FACTOR = "scaleFactor"

const BAYESIAN_BRIDGE_GIBBS_OP = "bayesianBridgeGibbsOperator"

const FACTOR_VALIDATION = "factorValidation"
const BIAS = "bias"
const SQUARED_ERROR = "squaredError"

const TRAIT_DATA_PROVIDER = "traitDataProvider"

const MULTIPLICATIVE_GAMMA_GIBBS = "multiplicativeGammaGibbsProvider"
const ROW_MULTIPLIERS = "rowMultipliers"
const ROW = "row"

const CACHED_MATRIX_INVERSE = "cachedMatrixInverse"
const COMPOUND_SYMMETRIC_MATRIX = "compoundSymmetricMatrix"
const AS_CORRELATION = "asCorrelation"
const IS_CHOLESKY = "isCholesky"
const DIAGONAL = "diagonal"
const OFF_DIAGONAL = "offDiagonal"

const LJK_CORRELATION_PRIOR = "LKJCorrelationPrior"
const SHAPE_PARAMETER = "shapeParameter"
const CHOLESKY = "cholesky"

const LOG_NORMAL_PRIOR = "logNormalPrior"
const MEAN_IN_REAL_SPACE = "meanInRealSpace"
const OFFSET = "offset"

const NORMAL_ORTHOGONAL_SUBSPACE = "normalOrthogonalSubspaceDistribution"

const NORMAL_MATRIX_NORM = "normalMatrixNormLikelihood"
const GLOBAL_PRECISION = "globalPrecision"
const MATRIX = "matrix"

const SCALED_MATRIX = "scaledMatrixParameter"
const MULTIVARIATE_GAMMA_LIKELIHOOD = "multivariateGammaLikelihood"

const MULTIPLICATIVE_PARAMETER = "multiplicativeParameter"

const GEODESIC_HMC = "geodesicHamiltonianMonteCarloOperator"
const NORMALIZED_LOADINGS_GRADIENT = "normalizedIntegratedLoadingsGradient"
const SCALE_LOADINGS_GRADIENT = "scaleIntegratedLoadingsGradient"

const MASK = "mask"
const MASKED_PARAMETER = "maskedParameter"

const INDICES = "indices"

const JOINT_OPERATOR = "jointOperator"
const INTEGRATED_FACTOR_GIBBS = "integratedFactorsGibbsOperator"
const SCALED_MATRIX_GRADIENT = "scaledMatrixGradient"
const COMPONENT = "component"
const SAMPLED_LOADINGS_GRADIENT = "sampledLoadingsGradient"

const PRECONDITIONING = "preconditioning"
const PRECONDITIONING_UPDATE_FREQUENCY = "preconditioningUpdateFrequency"
const PRECONDITIONING_MAX_UPDATE = "preconditioningMaxUpdate"
const PRECONDITIONING_DELAY = "preconditioningDelay"
const PRECONDITIONING_MEMORY = "preconditioningMemory"

const ADAPTIVE_DIAGONAL = "adaptiveDiagonal"
const MISSING_INDICATOR = "missingIndicator"

const LOCADINGS_SCALE_GIBBS = "loadingsScaleGibbsOperator"
const INDEPENDENT_NORMAL_DISTRIBUTION = "independentNormalDistributionModel"

const REJECTION_OPERATOR = "rejectionOperator"
const CONDITION = "condition"
const DESCENDING_AND_SPACED = "descendingAndSpaced"
const SPACING = "spacing"

const GENERAL_DATA_TYPE = "generalDataType"
const STATE = "state"
const CODE = "code"

const ALIGNMENT = "alignment"
const SEQUENCE = "sequence"
const PATTERNS = "patterns"

const FROM = "from"
const UNIQUE = "unique"

const ORDERED_LATENT_LIABILITY = "orderedLatentLiabilityLikelihood"
const TIP_TRAIT = "tipTrait"
const THRESHOLD = "threshold"
const NUM_CLASSES = "numClasses"

const COMPOUND_PARAMETER = "compoundParameter"

const EXTENDED_GIBBS_OPERATOR = "extendedLatentLiabilityGibbsOperator"

const N_TRAITS = "NTraits"
const N_DATA = "NData"

const FACTOR_PROPORTION_STATISTIC = "factorProportionStatistic"

end


text:
#PcPcOnSq(k::Int,i::Int,j::Int) = pc_on_sq[k,Int((i+1)*(i)/2)+(j)] # bug
###PcPcOnSq(k::Int,i::Int,j::Int) = pc_on_sq[Int((i-1)*(i)/2)+(j-1)+1,k] # also have bugs
#########PcPcOnSq(k::Int,i::Int,j::Int) = pc_on_sq[div(((i-1)*i),2)+(j-1)+1,k]
# pc_on_sq[k][(i)*((i)+1)/2+(j)]
PcPcOnSq(k::Int,i::Int,j::Int) = pc_on_sq[(div(i*(i+1),2)+j)+1,k+1]

function showPcPcOnSqIndex(k::Int, i::Int, j::Int)
    println("["*"$k,$i,$j"*"] = "*"pc_on_sq[",(i-1)*(i)/2+(j-1)+1,",",k,"]")
end

function probeEHash(p::Board, gs::GameStatus, key::UInt64)
    #key = hash(p.square) ⊻ hash(p.WhitePiecesInHands) ⊻ hash(p.BlackPiecesInHands) ⊻ hash(p.nextMove)
    contents = get( gs.ett, key, Int64(0xdeadcafe))
    if contents == Int64(0xdeadcafe)
        return false, 0
    else
        return true, contents
    end
end

function storeEHash(p::Board, gs::GameStatus, key::UInt64, value::Int64)
    gs.ett[key] = value
end

function make_list( list0::Array{Int,1}, list1::Array{Int,1}, p::Board, gs::GameStatus)
    list2::Array{Int,1} = [0 for x=1:35]::Array{Int,1}
    nlist::Int = 15
    score::Int = 0
    sq_bk0::Int = p.kingposB
    sq_wk0::Int = p.kingposW
    sq_bk1::Int = 82 - p.kingposW
    sq_wk1::Int = 82 - p.kingposB

    zerosOffset = leading_zeros(MaskOfBoard)

    # for i = 1:14
    #     println("list0[",i,"] = ", list0[i])
    # end
    # for i = 1:14
    #     println("list1[",i,"] = ", list1[i])
    # end

    score += kkp[1+ kkp_hand_pawn  + p.WhitePiecesInHands[MJFU], sq_bk0,sq_wk0]
    # println("wFU:",score,"kpp[",1+ kkp_hand_pawn  + p.WhitePiecesInHands[MJFU],",",sq_bk0,",",sq_wk0,"]")
    score += kkp[1+ kkp_hand_lance + p.WhitePiecesInHands[MJKY], sq_bk0,sq_wk0]
    score += kkp[1+ kkp_hand_knight+ p.WhitePiecesInHands[MJKE], sq_bk0,sq_wk0]
    score += kkp[1+ kkp_hand_silver+ p.WhitePiecesInHands[MJGI], sq_bk0,sq_wk0]
    score += kkp[1+ kkp_hand_gold  + p.WhitePiecesInHands[MJKI], sq_bk0,sq_wk0]
    score += kkp[1+ kkp_hand_bishop+ p.WhitePiecesInHands[MJKA], sq_bk0,sq_wk0]
    score += kkp[1+ kkp_hand_rook  + p.WhitePiecesInHands[MJHI], sq_bk0,sq_wk0]
    # println("wall:",score)
    score -= kkp[1+ kkp_hand_pawn  + p.BlackPiecesInHands[MJFU], sq_bk1,sq_wk1]
    #println("bFU:",score,"kpp[",1+ kkp_hand_pawn  + p.BlackPiecesInHands[MJFU],",",sq_bk1,",",sq_wk1,"]")
    score -= kkp[1+ kkp_hand_lance + p.BlackPiecesInHands[MJKY], sq_bk1,sq_wk1]
    score -= kkp[1+ kkp_hand_knight+ p.BlackPiecesInHands[MJKE], sq_bk1,sq_wk1]
    score -= kkp[1+ kkp_hand_silver+ p.BlackPiecesInHands[MJGI], sq_bk1,sq_wk1]
    score -= kkp[1+ kkp_hand_gold  + p.BlackPiecesInHands[MJKI], sq_bk1,sq_wk1]
    score -= kkp[1+ kkp_hand_bishop+ p.BlackPiecesInHands[MJKA], sq_bk1,sq_wk1]
    score -= kkp[1+ kkp_hand_rook  + p.BlackPiecesInHands[MJHI], sq_bk1,sq_wk1]

    # println("after add/sub kkp: ",score,"\n")

    n2::Int = 1
    bb::BitBoard = p.bb[MJFU]
    #DisplayBitBoard(bb,false)
    while bb > UInt128(0)
        sq::Int = trailing_zeros(bb)+1
        # println("sq=$sq, Inv(sq)=",(80 - sq+2))
        bb ⊻= BitSet[sq]
        list0[nlist] = f_pawn + sq - 1
        list2[n2]    = e_pawn + ((80 - sq)+2) - 1
        #println("list0[",nlist,"]=",f_pawn + sq + 1)
        #println("list2[",n2,"]=",e_pawn + ((80 - sq)+2))
        score += kkp[kkp_pawn + sq,sq_bk0,sq_wk0]
        #println("kkp[",(kkp_pawn+sq),"][",sq_bk0,"][",sq_wk0,"] = ", kkp[kkp_pawn + sq,sq_bk0,sq_wk0])
        nlist += 1
        n2    += 1
    end

    # println("after MJFU: score= ",score)

    bb = p.bb[MJGOFU]
    while bb > UInt128(0)
        sq::Int = trailing_zeros(bb)+1
        # println("sq=$sq, Inv(sq)=",(80 - sq+2))
        bb ⊻= BitSet[sq]
        list0[nlist] = e_pawn + sq - 1
        list2[n2]    = f_pawn + ((80 - sq)+2) - 1
        # println("list0[",nlist,"]=",e_pawn + sq)
        # println("list2[",n2,"]=",f_pawn + ((80 - sq)+2))
        score -= kkp[kkp_pawn + (80 - sq)+2,sq_bk1,sq_wk1]
        # println("kkp[",(kkp_pawn+(80 - sq)+2),"][",sq_bk1,"][",sq_wk1,"] = ", -kkp[kkp_pawn + (80 - sq)+2,sq_bk1,sq_wk1])
        nlist += 1
        n2    += 1
    end

    for i = 1:(n2-1)
        list1[nlist-i] = list2[i]
        #println("list1[",nlist,"-",i," = ", nlist-i, "] = list2[",i,"]")
    end

    # println("after MJGOFU: score= ",score)

    n2 = 1
    bb = p.bb[MJKY]
    while bb > UInt128(0)
        sq::Int = trailing_zeros(bb)+1
        # println("sq=$sq, Inv(sq)=",(80 - sq+2))
        bb ⊻= BitSet[sq]
        list0[nlist] = f_lance + sq - 1
        list2[n2]    = e_lance + ((80 - sq)+2) - 1
        # println("list0[",nlist,"]=",f_lance + sq)
        # println("list2[",n2,"]=",e_lance + ((80 - sq)+2))
        score += kkp[kkp_lance + sq,sq_bk0,sq_wk0]
        # println("kkp[",(kkp_lance+sq),"][",sq_bk0,"][",sq_wk0,"] = ", kkp[kkp_lance + sq,sq_bk0,sq_wk0])
        nlist += 1
        n2    += 1
    end
    # println("after MJKY: score= ",score)
    bb = p.bb[MJGOKY]
    while bb > UInt128(0)
        sq::Int = trailing_zeros(bb)+1
        # println("sq=$sq, Inv(sq)=",(80 - sq+2))
        bb ⊻= BitSet[sq]
        list0[nlist] = e_lance + sq - 1
        list2[n2]    = f_lance + ((80 - sq)+2) - 1
        score -= kkp[kkp_lance + (80 - sq)+2,sq_bk1,sq_wk1]
        nlist += 1
        n2    += 1
    end

    for i = 1:(n2-1)
        list1[nlist-i] = list2[i]
        #println("list1[",nlist,"-",i," = ", nlist-i, "] = list2[",i,"]")
    end

    # println("after MJGOKY: score= ",score)

    n2 = 1
    bb = p.bb[MJKE]
    while bb > UInt128(0)
        sq::Int = trailing_zeros(bb)+1
        # println("sq=$sq, Inv(sq)=",(80 - sq+2))
        bb ⊻= BitSet[sq]
        list0[nlist] = f_knight + sq - 1
        list2[n2]    = e_knight + ((80 - sq)+2) - 1
        score += kkp[kkp_knight + sq,sq_bk0,sq_wk0]
        nlist += 1
        n2    += 1
    end

    # println("after MJKE: score= ",score)

    bb = p.bb[MJGOKE]
    while bb > UInt128(0)
        sq::Int = trailing_zeros(bb)+1
        # println("sq=$sq, Inv(sq)=",(80 - sq+2))
        bb ⊻= BitSet[sq]
        list0[nlist] = e_knight + sq - 1
        list2[n2]    = f_knight + ((80 - sq)+2) - 1
        score -= kkp[kkp_knight + (80 - sq)+2,sq_bk1,sq_wk1]
        nlist += 1
        n2    += 1
    end

    for i = 1:(n2-1)
        list1[nlist-i] = list2[i]
        #println("list1[",nlist,"-",i," = ", nlist-i, "] = list2[",i,"]")
    end

    # println("after MJGOKE: score= ",score)

    n2 = 1
    bb = p.bb[MJGI]
    while bb > UInt128(0)
        sq::Int = trailing_zeros(bb)+1
        # println("sq=$sq, Inv(sq)=",(80 - sq+2))
        bb ⊻= BitSet[sq]
        list0[nlist] = f_silver + sq - 1
        list2[n2]    = e_silver + ((80 - sq)+2) - 1
        score += kkp[kkp_silver + sq,sq_bk0,sq_wk0]
        nlist += 1
        n2    += 1
    end

    # println("after MJGI: score= ",score)

    bb = p.bb[MJGOGI]
    while bb > UInt128(0)
        sq::Int = trailing_zeros(bb)+1
        # println("sq=$sq, Inv(sq)=",(80 - sq+2))
        bb ⊻= BitSet[sq]
        list0[nlist] = e_silver + sq - 1
        list2[n2]    = f_silver + ((80 - sq)+2) - 1
        score -= kkp[kkp_silver + (80 - sq)+2, sq_bk1,sq_wk1]
        nlist += 1
        n2    += 1
    end

    for i = 1:(n2-1)
        list1[nlist-i] = list2[i]
        #println("list1[",nlist,"-",i," = ", nlist-i, "] = list2[",i,"]")
    end
    # println("after MJGOGI: score= ",score)
    n2 = 1
    bb = (p.bb[MJKI]|p.bb[MJTO]|p.bb[MJNY]|p.bb[MJNK]|p.bb[MJNG])
    while bb > UInt128(0)
        sq::Int = trailing_zeros(bb)+1
        # println("sq=$sq, Inv(sq)=",(80 - sq+2))
        bb ⊻= BitSet[sq]
        list0[nlist] = f_gold + sq - 1
        list2[n2]    = e_gold + ((80 - sq)+2) - 1
        score += kkp[kkp_gold + sq,sq_bk0,sq_wk0]
        nlist += 1
        n2    += 1
    end
    # println("after MJKI etc: score= ",score)
    bb = (p.bb[MJGOKI]|p.bb[MJGOTO]|p.bb[MJGONY]|p.bb[MJGONK]|p.bb[MJGONG])
    while bb > UInt128(0)
        sq::Int = trailing_zeros(bb)+1
        # println("sq=$sq, Inv(sq)=",(80 - sq+2))
        bb ⊻= BitSet[sq]
        list0[nlist] = e_gold + sq - 1
        list2[n2]    = f_gold + ((80 - sq)+2) - 1
        score -= kkp[kkp_gold + (80 - sq)+2,sq_bk1,sq_wk1]
        nlist += 1
        n2    += 1
    end

    for i = 1:(n2-1)
        list1[nlist-i] = list2[i]
        #println("list1[",nlist,"-",i," = ", nlist-i, "] = list2[",i,"]")
    end
    # println("after MJGOKI etc: score= ",score)
    n2 = 1
    bb = p.bb[MJKA]
    while bb > UInt128(0)
        sq::Int = trailing_zeros(bb)+1
        # println("sq=$sq, Inv(sq)=",(80 - sq+2))
        bb ⊻= BitSet[sq]
        list0[nlist] = f_bishop + sq - 1
        list2[n2]    = e_bishop + ((80 - sq)+2) - 1
        score += kkp[kkp_bishop + sq,sq_bk0,sq_wk0]
        nlist += 1
        n2    += 1
    end
    # println("after MJKA etc: score= ",score)
    bb = p.bb[MJGOKA]
    while bb > UInt128(0)
        sq::Int = trailing_zeros(bb)+1
        # println("sq=$sq, Inv(sq)=",(80 - sq+2))
        bb ⊻= BitSet[sq]
        list0[nlist] = e_bishop + sq - 1
        list2[n2]    = f_bishop + ((80 - sq)+2) - 1
        score -= kkp[kkp_bishop + (80 - sq)+2,sq_bk1,sq_wk1]
        nlist += 1
        n2    += 1
    end

    for i = 1:(n2-1)
        list1[nlist-i] = list2[i]
        #println("list1[",nlist,"-",i," = ", nlist-i, "] = list2[",i,"]")
    end
    # println("after MJGOKA etc: score= ",score)
    n2 = 1
    bb = p.bb[MJUM]
    while bb > UInt128(0)
        sq::Int = trailing_zeros(bb)+1
        # println("sq=$sq, Inv(sq)=",(80 - sq+2))
        bb ⊻= BitSet[sq]
        list0[nlist] = f_horse + sq - 1
        list2[n2]    = e_horse + ((80 - sq)+2) - 1
        score += kkp[kkp_horse + sq,sq_bk0,sq_wk0]
        nlist += 1
        n2    += 1
    end
    # println("after MJUM etc: score= ",score)

    bb = p.bb[MJGOUM]
    while bb > UInt128(0)
        sq::Int = trailing_zeros(bb)+1
        # println("sq=$sq, Inv(sq)=",(80 - sq+2))
        bb ⊻= BitSet[sq]
        list0[nlist] = e_horse + sq - 1
        list2[n2]    = f_horse + ((80 - sq)+2) - 1
        score -= kkp[kkp_horse + (80 - sq)+2,sq_bk1,sq_wk1]
        nlist += 1
        n2    += 1
    end

    for i = 1:(n2-1)
        list1[nlist-i] = list2[i]
        #println("list1[",nlist,"-",i," = ", nlist-i, "] = list2[",i,"]")
    end
    # println("after MJGOUM etc: score= ",score)
    n2 = 1
    bb = p.bb[MJHI]
    while bb > UInt128(0)
        sq::Int = trailing_zeros(bb)+1
        # println("sq=$sq, Inv(sq)=",(80 - sq+2))
        bb ⊻= BitSet[sq]
        list0[nlist] = f_rook + sq - 1
        list2[n2]    = e_rook + ((80 - sq)+2) - 1
        score += kkp[kkp_rook + sq,sq_bk0,sq_wk0]
        nlist += 1
        n2    += 1
    end
    # println("after MJHI etc: score= ",score)
    bb = p.bb[MJGOHI]
    while bb > UInt128(0)
        sq::Int = trailing_zeros(bb)+1
        # println("sq=$sq, Inv(sq)=",(80 - sq+2))
        bb ⊻= BitSet[sq]
        list0[nlist] = e_rook + sq - 1
        list2[n2]    = f_rook + ((80 - sq)+2) - 1
        score -= kkp[kkp_rook + (80 - sq)+2,sq_bk1,sq_wk1]
        nlist += 1
        n2    += 1
    end

    for i = 1:(n2-1)
        list1[nlist-i] = list2[i]
        #println("list1[",nlist,"-",i," = ", nlist-i, "] = list2[",i,"]")
    end
    # println("after MJGOHI etc: score= ",score)
    n2 = 1
    bb = p.bb[MJRY]
    while bb > UInt128(0)
        sq::Int = trailing_zeros(bb)+1
        # println("sq=$sq, Inv(sq)=",(80 - sq+2))
        bb ⊻= BitSet[sq]
        list0[nlist] = f_dragon + sq - 1
        list2[n2]    = e_dragon + ((80 - sq)+2) - 1
        score += kkp[kkp_dragon + sq,sq_bk0,sq_wk0]
        nlist += 1
        n2    += 1
    end
    # println("after MJRY etc: score= ",score)
    bb = p.bb[MJGORY]
    while bb > UInt128(0)
        sq::Int = trailing_zeros(bb)+1
        # println("sq=$sq, Inv(sq)=",(80 - sq+2))
        bb ⊻= BitSet[sq]
        list0[nlist] = e_dragon + sq - 1
        list2[n2]    = f_dragon + ((80 - sq)+2) - 1
        score -= kkp[kkp_dragon + (80 - sq)+2,sq_bk1,sq_wk1]
        nlist += 1
        n2    += 1
    end

    for i = 1:(n2-1)
        list1[nlist-i] = list2[i]
        #println("list1[",nlist,"-",i," = ", nlist-i, "] = list2[",i,"]")
    end
    # println("after MJGORY etc: score= ",score)
    if nlist > 53
        println("nlist is larger than 53.")
    end
    #println("makelist score =", score)
    return nlist-1, score
end

function EvalBonanza(nextMove::Int, p::Board, gs::GameStatus)
    # list0::Array{Int,1} = [0 for x=1:53]::Array{Int,1}
    # list1::Array{Int,1} = [0 for x=1:53]::Array{Int,1}
    list0 = zeros(Int, 54)::Array{Int,1}
    list1 = zeros(Int, 54)::Array{Int,1}

    # hash lookup!

    key = hash(p.square) ⊻ hash(p.WhitePiecesInHands) ⊻ hash(p.BlackPiecesInHands) ⊻ hash(p.nextMove)

    inHash, evalue = probeEHash(p, gs, key)
    if inHash == true
        if nextMove == GOTE
            evalue = -evalue
        end

        evalue = Int64(div(evalue,32))

        noise = (rand(Int64) % 10) - 5
        evalue = evalue + noise

        return evalue
    end
    #println("list0[1] = ", f_hand_pawn, "+",p.WhitePiecesInHands[MJFU])
    list0[ 1] = f_hand_pawn   + p.WhitePiecesInHands[MJFU]
    #println("list0[2] = ", e_hand_pawn, "+",p.BlackPiecesInHands[MJFU])
    list0[ 2] = e_hand_pawn   + p.BlackPiecesInHands[MJFU]

    list0[ 3] = f_hand_lance  + p.WhitePiecesInHands[MJKY]
    list0[ 4] = e_hand_lance  + p.BlackPiecesInHands[MJKY]
    list0[ 5] = f_hand_knight + p.WhitePiecesInHands[MJKE]
    list0[ 6] = e_hand_knight + p.BlackPiecesInHands[MJKE]
    list0[ 7] = f_hand_silver + p.WhitePiecesInHands[MJGI]
    list0[ 8] = e_hand_silver + p.BlackPiecesInHands[MJGI]
    list0[ 9] = f_hand_gold   + p.WhitePiecesInHands[MJKI]
    list0[10] = e_hand_gold   + p.BlackPiecesInHands[MJKI]
    list0[11] = f_hand_bishop + p.WhitePiecesInHands[MJKA]
    list0[12] = e_hand_bishop + p.BlackPiecesInHands[MJKA]
    list0[13] = f_hand_rook   + p.WhitePiecesInHands[MJHI]
    list0[14] = e_hand_rook   + p.BlackPiecesInHands[MJHI]

    list1[ 1] = f_hand_pawn   + p.BlackPiecesInHands[MJFU]
    list1[ 2] = e_hand_pawn   + p.WhitePiecesInHands[MJFU]
    list1[ 3] = f_hand_lance  + p.BlackPiecesInHands[MJKY]
    list1[ 4] = e_hand_lance  + p.WhitePiecesInHands[MJKY]
    list1[ 5] = f_hand_knight + p.BlackPiecesInHands[MJKE]
    list1[ 6] = e_hand_knight + p.WhitePiecesInHands[MJKE]
    list1[ 7] = f_hand_silver + p.BlackPiecesInHands[MJGI]
    list1[ 8] = e_hand_silver + p.WhitePiecesInHands[MJGI]
    list1[ 9] = f_hand_gold   + p.BlackPiecesInHands[MJKI]
    list1[10] = e_hand_gold   + p.WhitePiecesInHands[MJKI]
    list1[11] = f_hand_bishop + p.BlackPiecesInHands[MJKA]
    list1[12] = e_hand_bishop + p.WhitePiecesInHands[MJKA]
    list1[13] = f_hand_rook   + p.BlackPiecesInHands[MJHI]
    list1[14] = e_hand_rook   + p.WhitePiecesInHands[MJHI]

    nlist::Int = 0
    score::Int = 0

    nlist, score = make_list( list0, list1, p, gs)
    #for i = 1:nlist
    #    println("list0[",i-1,"] = ", list0[i])
    #end
    #for i = 1:nlist
    #    println("list1[",i-1,"] = ", list1[i])
    #end

    #for i = 1:nlist
    #    if list0[i] >= fe_end
    #        println("list0[",i,"]=",list0[i])
    #    end
    #    if list1[i] >= fe_end
    #        println("list1[",i,"]=",list1[i])
    #    end
    #end
    sq_bk::Int = p.kingposW - 1
    sq_wk::Int = 82 - p.kingposB - 1
    #println("sq_bk=$(sq_bk), sq_wk = $(sq_wk)")
    #println("p.kingposW = $(p.kingposW), p.kingposB = $(p.kingposB)")
    sum::Int = 0

    for i = 0:(nlist-1)
        k0::Int = list0[i+1]
        k1::Int = list1[i+1]

        for j = 0:(i)
            #println("(i,j)=($i,$j)")
            l0::Int = list0[j+1]
            l1::Int = list1[j+1]
            #try
            sum = sum + PcPcOnSq( sq_bk, k0, l0)
            #showPcPcOnSqIndex(sq_bk, k0, l0)
            #println("sq_bk=",sq_bk,",k0=",k0,",l0=",l0)

            sum = sum - PcPcOnSq( sq_wk, k1, l1)
            #showPcPcOnSqIndex(sq_wk, k1, l1)
            #println("sq_wk=",sq_wk,",k1=",k1,",l1=",l1)
            #catch
            #    quit()
            #    println("sq_bk=",sq_bk,",sq_wk=",sq_wk)
            #    println("k0=",k0,",k1=",k1)
            #    println("l0=",l0,",l1=",l1)
            #    println("i,j=",i,",",j)
            #    quit()
            #end
        end
    end

    score = score + sum
    eva = Eval( SENTE, p, gs)
    #println("MATERIAL(J) = ", eva)
    #println("sum(J) = ", sum)
    score = score + (32 * eva)

    # store eHash!
    if inHash == false # always false
        storeEHash(p, gs, key, Int64(score))
    end

    # println("Eval=", eva, ", Score=", score)

    if nextMove == GOTE
        score = -score
    end

    score = Int(div(score,32))

    #noise = (rand(Int64) % 10) - 5
    #score += noise

    #println("XXX")
    #DisplayBoard(p)
    #println("score = $(score)")

    return score
end

ext: jl

max_stars_repo_name: kimrin/MechajyoForWindows

max_issues_repo_name: kimrin/MechajyoForWindows

max_forks_repo_name: kimrin/MechajyoForWindows


text:
# FIXME: handle extern
const randombytes_internal_implementation = randombytes_implementation
const randombytes_salsa20_implementation = randombytes_internal_implementation


text:
promote_rule(::Type{MarkInt128}, ::Type{MarkInt64}) = MarkInt128
promote_rule(::Type{MarkInt128}, ::Type{MarkInt32}) = MarkInt128
promote_rule(::Type{MarkInt128}, ::Type{MarkInt16}) = MarkInt128
promote_rule(::Type{MarkInt128}, ::Type{MarkInt8}) = MarkInt128
promote_rule(::Type{MarkInt64}, ::Type{MarkInt32}) = MarkInt64
promote_rule(::Type{MarkInt64}, ::Type{MarkInt16}) = MarkInt64
promote_rule(::Type{MarkInt64}, ::Type{MarkInt8}) = MarkInt64
promote_rule(::Type{MarkInt32}, ::Type{MarkInt16}) = MarkInt32
promote_rule(::Type{MarkInt32}, ::Type{MarkInt8}) = MarkInt32
promote_rule(::Type{MarkInt16}, ::Type{MarkInt8}) = MarkInt16

convert(::Type{MarkInt128}, x::MarkInt64) = isunmarked(x) ? reinterpret(MarkInt128, Int128(ityped(x)>>1)<<1) : mark(reinterpret(MarkInt128, Int128(ityped(x)>>1)<<1))
convert(::Type{MarkInt128}, x::MarkInt32) = isunmarked(x) ? reinterpret(MarkInt128, Int128(ityped(x)>>1)<<1) : mark(reinterpret(MarkInt128, Int128(ityped(x)>>1)<<1))
convert(::Type{MarkInt128}, x::MarkInt16) = isunmarked(x) ? reinterpret(MarkInt128, Int128(ityped(x)>>1)<<1) : mark(reinterpret(MarkInt128, Int128(ityped(x)>>1)<<1))
convert(::Type{MarkInt128}, x::MarkInt8) = isunmarked(x) ? reinterpret(MarkInt128, Int128(ityped(x)>>1)<<1) : mark(reinterpret(MarkInt128, Int128(ityped(x)>>1)<<1))
convert(::Type{MarkInt64}, x::MarkInt32) = isunmarked(x) ? reinterpret(MarkInt64, Int64(ityped(x)>>1)<<1) : mark(reinterpret(MarkInt64, Int64(ityped(x)>>1)<<1))
convert(::Type{MarkInt64}, x::MarkInt16) = isunmarked(x) ? reinterpret(MarkInt64, Int64(ityped(x)>>1)<<1) : mark(reinterpret(MarkInt64, Int64(ityped(x)>>1)<<1))
convert(::Type{MarkInt64}, x::MarkInt8) = isunmarked(x) ? reinterpret(MarkInt64, Int64(ityped(x)>>1)<<1) : mark(reinterpret(MarkInt64, Int64(ityped(x)>>1)<<1))
convert(::Type{MarkInt32}, x::MarkInt16) = isunmarked(x) ? reinterpret(MarkInt32, Int32(ityped(x)>>1)<<1) : mark(reinterpret(MarkInt32, Int32(ityped(x)>>1)<<1))
convert(::Type{MarkInt32}, x::MarkInt8) = isunmarked(x) ? reinterpret(MarkInt32, Int32(ityped(x)>>1)<<1) : mark(reinterpret(MarkInt32, Int32(ityped(x)>>1)<<1))
convert(::Type{MarkInt16}, x::MarkInt8) = isunmarked(x) ? reinterpret(MarkInt16, Int16(ityped(x)>>1)<<1) : mark(reinterpret(MarkInt16, Int16(ityped(x)>>1)<<1))

promote_rule(::Type{MarkUInt128}, ::Type{MarkUInt64}) = MarkUInt128
promote_rule(::Type{MarkUInt128}, ::Type{MarkUInt32}) = MarkUInt128
promote_rule(::Type{MarkUInt128}, ::Type{MarkUInt16}) = MarkUInt128
promote_rule(::Type{MarkUInt128}, ::Type{MarkUInt8}) = MarkUInt128
promote_rule(::Type{MarkUInt64}, ::Type{MarkUInt32}) = MarkUInt64
promote_rule(::Type{MarkUInt64}, ::Type{MarkUInt16}) = MarkUInt64
promote_rule(::Type{MarkUInt64}, ::Type{MarkUInt8}) = MarkUInt64
promote_rule(::Type{MarkUInt32}, ::Type{MarkUInt16}) = MarkUInt32
promote_rule(::Type{MarkUInt32}, ::Type{MarkUInt8}) = MarkUInt32
promote_rule(::Type{MarkUInt16}, ::Type{MarkUInt8}) = MarkUInt16

convert(::Type{MarkUInt128}, x::MarkUInt64) = isunmarked(x) ? reinterpret(MarkUInt128, UInt128(ityped(x)>>1)<<1) : mark(reinterpret(MarkUInt128, UInt128(ityped(x)>>1)<<1))
convert(::Type{MarkUInt128}, x::MarkUInt32) = isunmarked(x) ? reinterpret(MarkUInt128, UInt128(ityped(x)>>1)<<1) : mark(reinterpret(MarkUInt128, UInt128(ityped(x)>>1)<<1))
convert(::Type{MarkUInt128}, x::MarkUInt16) = isunmarked(x) ? reinterpret(MarkUInt128, UInt128(ityped(x)>>1)<<1) : mark(reinterpret(MarkUInt128, UInt128(ityped(x)>>1)<<1))
convert(::Type{MarkUInt128}, x::MarkUInt8) = isunmarked(x) ? reinterpret(MarkUInt128, UInt128(ityped(x)>>1)<<1) : mark(reinterpret(MarkUInt128, UInt128(ityped(x)>>1)<<1))
convert(::Type{MarkUInt64}, x::MarkUInt32) = isunmarked(x) ? reinterpret(MarkUInt64, UInt64(ityped(x)>>1)<<1) : mark(reinterpret(MarkUInt64, UInt64(ityped(x)>>1)<<1))
convert(::Type{MarkUInt64}, x::MarkUInt16) = isunmarked(x) ? reinterpret(MarkUInt64, UInt64(ityped(x)>>1)<<1) : mark(reinterpret(MarkUInt64, UInt64(ityped(x)>>1)<<1))
convert(::Type{MarkUInt64}, x::MarkUInt8) = isunmarked(x) ? reinterpret(MarkUInt64, UInt64(ityped(x)>>1)<<1) : mark(reinterpret(MarkUInt64, UInt64(ityped(x)>>1)<<1))
convert(::Type{MarkUInt32}, x::MarkUInt16) = isunmarked(x) ? reinterpret(MarkUInt32, UInt32(ityped(x)>>1)<<1) : mark(reinterpret(MarkUInt32, UInt32(ityped(x)>>1)<<1))
convert(::Type{MarkUInt32}, x::MarkUInt8) = isunmarked(x) ? reinterpret(MarkUInt32, UInt32(ityped(x)>>1)<<1) : mark(reinterpret(MarkUInt32, UInt32(ityped(x)>>1)<<1))
convert(::Type{MarkUInt16}, x::MarkUInt8) = isunmarked(x) ? reinterpret(MarkUInt16, UInt16(ityped(x)>>1)<<1) : mark(reinterpret(MarkUInt16, UInt16(ityped(x)>>1)<<1))

promote_rule(::Type{MarkInt128}, ::Type{Int128}) = MarkInt128
promote_rule(::Type{MarkInt128}, ::Type{Int64}) = MarkInt128
promote_rule(::Type{MarkInt128}, ::Type{Int32}) = MarkInt128
promote_rule(::Type{MarkInt128}, ::Type{Int16}) = MarkInt128
promote_rule(::Type{MarkInt128}, ::Type{Int8}) = MarkInt128
promote_rule(::Type{MarkInt64}, ::Type{Int64}) = MarkInt64
promote_rule(::Type{MarkInt64}, ::Type{Int32}) = MarkInt64
promote_rule(::Type{MarkInt64}, ::Type{Int16}) = MarkInt64
promote_rule(::Type{MarkInt64}, ::Type{Int8}) = MarkInt64
promote_rule(::Type{MarkInt32}, ::Type{Int32}) = MarkInt32
promote_rule(::Type{MarkInt32}, ::Type{Int16}) = MarkInt32
promote_rule(::Type{MarkInt32}, ::Type{Int8}) = MarkInt32
promote_rule(::Type{MarkInt16}, ::Type{Int16}) = MarkInt16
promote_rule(::Type{MarkInt16}, ::Type{Int8}) = MarkInt16
promote_rule(::Type{MarkInt8}, ::Type{Int8}) = MarkInt8

convert(::Type{MarkInt128}, x::Int128) = Unmarked(x)
convert(::Type{MarkInt128}, x::Int64) = Unmarked(Int128(x))
convert(::Type{MarkInt128}, x::Int32) = Unmarked(Int128(x))
convert(::Type{MarkInt128}, x::Int16) = Unmarked(Int128(x))
convert(::Type{MarkInt128}, x::Int8) = Unmarked(Int128(x))
convert(::Type{MarkInt64}, x::Int64) = Unmarked(x)
convert(::Type{MarkInt64}, x::Int32) = Unmarked(Int64(x))
convert(::Type{MarkInt64}, x::Int16) = Unmarked(Int64(x))
convert(::Type{MarkInt64}, x::Int8) = Unmarked(Int64(x))
convert(::Type{MarkInt32}, x::Int32) = Unmarked(x)
convert(::Type{MarkInt32}, x::Int16) = Unmarked(Int32(x))
convert(::Type{MarkInt32}, x::Int8) = Unmarked(Int32(x))
convert(::Type{MarkInt16}, x::Int16) = Unmarked(x)
convert(::Type{MarkInt16}, x::Int8) = Unmarked(Int16(x))
convert(::Type{MarkInt8}, x::Int8) = Unmarked(x)

promote_rule(::Type{MarkUInt128}, ::Type{UInt128}) = MarkUInt128
promote_rule(::Type{MarkUInt128}, ::Type{UInt64}) = MarkUInt128
promote_rule(::Type{MarkUInt128}, ::Type{UInt32}) = MarkUInt128
promote_rule(::Type{MarkUInt128}, ::Type{UInt16}) = MarkUInt128
promote_rule(::Type{MarkUInt128}, ::Type{UInt8}) = MarkUInt128
promote_rule(::Type{MarkUInt64}, ::Type{UInt64}) = MarkUInt64
promote_rule(::Type{MarkUInt64}, ::Type{UInt32}) = MarkUInt64
promote_rule(::Type{MarkUInt64}, ::Type{UInt16}) = MarkUInt64
promote_rule(::Type{MarkUInt64}, ::Type{UInt8}) = MarkUInt64
promote_rule(::Type{MarkUInt32}, ::Type{UInt32}) = MarkUInt32
promote_rule(::Type{MarkUInt32}, ::Type{UInt16}) = MarkUInt32
promote_rule(::Type{MarkUInt32}, ::Type{UInt8}) = MarkUInt32
promote_rule(::Type{MarkUInt16}, ::Type{UInt16}) = MarkUInt16
promote_rule(::Type{MarkUInt16}, ::Type{UInt8}) = MarkUInt16
promote_rule(::Type{MarkUInt8}, ::Type{UInt8}) = MarkUInt8

convert(::Type{MarkUInt128}, x::UInt128) = Unmarked(x)
convert(::Type{MarkUInt128}, x::UInt64) = Unmarked(UInt128(x))
convert(::Type{MarkUInt128}, x::UInt32) = Unmarked(UInt128(x))
convert(::Type{MarkUInt128}, x::UInt16) = Unmarked(UInt128(x))
convert(::Type{MarkUInt128}, x::UInt8) = Unmarked(UInt128(x))
convert(::Type{MarkUInt64}, x::UInt64) = Unmarked(x)
convert(::Type{MarkUInt64}, x::UInt32) = Unmarked(UInt64(x))
convert(::Type{MarkUInt64}, x::UInt16) = Unmarked(UInt64(x))
convert(::Type{MarkUInt64}, x::UInt8) = Unmarked(UInt64(x))
convert(::Type{MarkUInt32}, x::UInt32) = Unmarked(x)
convert(::Type{MarkUInt32}, x::UInt16) = Unmarked(UInt32(x))
convert(::Type{MarkUInt32}, x::UInt8) = Unmarked(UInt32(x))
convert(::Type{MarkUInt16}, x::UInt16) = Unmarked(x)
convert(::Type{MarkUInt16}, x::UInt8) = Unmarked(UInt16(x))
convert(::Type{MarkUInt8}, x::UInt8) = Unmarked(x)

MarkInt128(x::MarkInt64) = ismarked(x) ? Marked(Int128(itype(x))) : Unmarked(Int128(itype(x)))
MarkInt128(x::MarkInt32) = ismarked(x) ? Marked(Int128(itype(x))) : Unmarked(Int128(itype(x)))
MarkInt128(x::MarkInt16) = ismarked(x) ? Marked(Int128(itype(x))) : Unmarked(Int128(itype(x)))
MarkInt128(x::MarkInt8) = ismarked(x) ? Marked(Int128(itype(x))) : Unmarked(Int128(itype(x)))
MarkInt64(x::MarkInt32) = ismarked(x) ? Marked(Int64(itype(x))) : Unmarked(Int64(itype(x)))
MarkInt64(x::MarkInt16) = ismarked(x) ? Marked(Int64(itype(x))) : Unmarked(Int64(itype(x)))
MarkInt64(x::MarkInt8) = ismarked(x) ? Marked(Int64(itype(x))) : Unmarked(Int64(itype(x)))
MarkInt32(x::MarkInt16) = ismarked(x) ? Marked(Int32(itype(x))) : Unmarked(Int32(itype(x)))
MarkInt32(x::MarkInt8) = ismarked(x) ? Marked(Int32(itype(x))) : Unmarked(Int32(itype(x)))
MarkInt16(x::MarkInt8) = ismarked(x) ? Marked(Int16(itype(x))) : Unmarked(Int16(itype(x)))

MarkUInt128(x::MarkUInt64) = ismarked(x) ? Marked(UInt128(itype(x))) : Unmarked(UInt128(itype(x)))
MarkUInt128(x::MarkUInt32) = ismarked(x) ? Marked(UInt128(itype(x))) : Unmarked(UInt128(itype(x)))
MarkUInt128(x::MarkUInt16) = ismarked(x) ? Marked(UInt128(itype(x))) : Unmarked(UInt128(itype(x)))
MarkUInt128(x::MarkUInt8) = ismarked(x) ? Marked(UInt128(itype(x))) : Unmarked(UInt128(itype(x)))
MarkUInt64(x::MarkUInt32) = ismarked(x) ? Marked(UInt64(itype(x))) : Unmarked(UInt64(itype(x)))
MarkUInt64(x::MarkUInt16) = ismarked(x) ? Marked(UInt64(itype(x))) : Unmarked(UInt64(itype(x)))
MarkUInt64(x::MarkUInt8) = ismarked(x) ? Marked(UInt64(itype(x))) : Unmarked(UInt64(itype(x)))
MarkUInt32(x::MarkUInt16) = ismarked(x) ? Marked(UInt32(itype(x))) : Unmarked(UInt32(itype(x)))
MarkUInt32(x::MarkUInt8) = ismarked(x) ? Marked(UInt32(itype(x))) : Unmarked(UInt32(itype(x)))
MarkUInt16(x::MarkUInt8) = ismarked(x) ? Marked(UInt16(itype(x))) : Unmarked(UInt16(itype(x)))

MarkInt128(x::I) where {I<:Union{Int64, Int32, Int16, Int8}} = MarkInt128(Int128(x))
MarkUInt128(x::I) where {I<:Union{UInt64, UInt32, UInt16, UInt8}} = MarkUInt128(UInt128(x))
MarkInt64(x::I) where {I<:Union{Int32, Int16, Int8}} = MarkInt64(Int64(x))
MarkUInt64(x::I) where {I<:Union{UInt32, UInt16, UInt8}} = MarkUInt64(UInt64(x))
MarkInt32(x::I) where {I<:Union{Int16, Int8}} = MarkInt32(Int32(x))
MarkUInt32(x::I) where {I<:Union{UInt16, UInt8}} = MarkUInt32(UInt32(x))
MarkInt16(x::Int8) = MarkInt16(Int16(x))
MarkUInt16(x::UInt8) = MarkUInt16(UInt16(x))

convert(::Type{Int64}, x::Int128) = Int64(x)
convert(::Type{Int64}, x::Int64) = x
convert(::Type{Int64}, x::Int32) = Int64(x)
convert(::Type{Int64}, x::Int16) = Int64(x)
convert(::Type{Int64}, x::Int8) = Int64(x)
convert(::Type{Int32}, x::Int128) = Int32(x)
convert(::Type{Int32}, x::Int64) = Int32(x)
convert(::Type{Int32}, x::Int32) = x
convert(::Type{Int32}, x::Int16) = Int32(x)
convert(::Type{Int32}, x::Int8) = Int32(x)
convert(::Type{Int16}, x::Int128) = Int16(x)
convert(::Type{Int16}, x::Int64) = Int16(x)
convert(::Type{Int16}, x::Int32) = Int16(x)
convert(::Type{Int16}, x::Int16) = x
convert(::Type{Int16}, x::Int8) = Int16(x)

convert(::Type{UInt64}, x::UInt128) = UInt64(x)
convert(::Type{UInt64}, x::UInt64) = x
convert(::Type{UInt64}, x::UInt32) = UInt64(x)
convert(::Type{UInt64}, x::UInt16) = UInt64(x)
convert(::Type{UInt64}, x::UInt8) = UInt64(x)
convert(::Type{UInt32}, x::UInt128) = UInt32(x)
convert(::Type{UInt32}, x::UInt64) = UInt32(x)
convert(::Type{UInt32}, x::UInt32) = x
convert(::Type{UInt32}, x::UInt16) = UInt32(x)
convert(::Type{UInt32}, x::UInt8) = UInt32(x)
convert(::Type{UInt16}, x::UInt128) = UInt16(x)
convert(::Type{UInt16}, x::UInt64) = UInt16(x)
convert(::Type{UInt16}, x::UInt32) = UInt16(x)
convert(::Type{UInt16}, x::UInt16) = x
convert(::Type{UInt16}, x::UInt8) = UInt16(x)


text:
function test_header(hdStr :: String)
    test_divider();
    println("Testing  $hdStr")
    return nothing
end

function test_divider()
    println("-------------------")
    return nothing
end

# Put files generated by tests here
function test_dir()
	return "~/Documents/julia/test_files"
end

# ----------


text:
# using Revise, Plots
using Test
	using BifurcationKit, Parameters, Setfield, LinearAlgebra, ForwardDiff, SparseArrays
	const BK = BifurcationKit
##################################################################
# The goal of these tests is to test all combinaisons of options
##################################################################

norminf = x -> norm(x, Inf)

function Fsl!(f, u, p, t)
	@unpack r, μ, ν, c3 = p
	u1 = u[1]
	u2 = u[2]

	ua = u1^2 + u2^2

	f[1] = r * u1 - ν * u2 - ua * (c3 * u1 - μ * u2)
	f[2] = r * u2 + ν * u1 - ua * (c3 * u2 + μ * u1)
	return f
end

Fsl(x, p) = Fsl!(similar(x), x, p, 0.)

####################################################################################################
par_sl = (r = 0.5, μ = 0., ν = 1.0, c3 = 1.0)
u0 = [.001, .001]
par_hopf = (@set par_sl.r = 0.1)
####################################################################################################
# continuation, Hopf bifurcation point detection
optconteq = ContinuationPar(ds = -0.01, detectBifurcation = 3, pMin = -0.5, nInversion = 8)
br, = continuation(Fsl, u0, par_hopf, (@lens _.r), optconteq)
####################################################################################################
poTrap = PeriodicOrbitTrapProblem(
	Fsl, (x, p) -> sparse(ForwardDiff.jacobian(z -> Fsl(z, p), x)), # we put sparse to try the different linear solvers
	[1., 0.],
	zeros(2),
	10, 2)

BK.hasHessian(poTrap)
show(poTrap)

# guess for the periodic orbit
orbitguess_f = reduce(vcat, [√(par_hopf.r) .* [cos(θ), sin(θ)] for θ in LinRange(0, 2pi, poTrap.M)])
	push!(orbitguess_f, 2pi)
optn_po = NewtonPar()
opts_po_cont = ContinuationPar(dsmax = 0.02, ds = 0.001, pMax = 2.2, maxSteps = 25, newtonOptions = optn_po, saveSolEveryStep = 1)

lsdef = DefaultLS()
lsit = GMRESKrylovKit()
for (ind, jacobianPO) in enumerate([:Dense, :FullLU, :BorderedLU, :FullSparseInplace, :BorderedSparseInplace, :FullMatrixFree, :BorderedMatrixFree])
	@show jacobianPO, ind
	_ls = ind > 5 ? lsit : lsdef
	outpo_f, _, flag = newton(poTrap,
		orbitguess_f, par_hopf, (@set optn_po.linsolver = _ls);
		jacobianPO = jacobianPO,
		normN = norminf)
	@test flag

	br_po, = continuation(poTrap, outpo_f,
		par_hopf, (@lens _.r),	(@set opts_po_cont.newtonOptions.linsolver = _ls), jacobianPO = jacobianPO;
		verbosity = 0,	plot = false,
		# plotSolution = (x, p; kwargs...) -> BK.plotPeriodicPOTrap(x, poTrap.M, 2, 1; ratio = 2, kwargs...),
		printSolution = (u, p) -> BK.getAmplitude(poTrap, u, par_hopf; ratio = 1), normC = norminf)

	BK.getPeriodicOrbit(br_po, 1)
end

outpo_f, = newton(poTrap, orbitguess_f, par_hopf, optn_po; jacobianPO = :Dense)
outpo = reshape(outpo_f[1:end-1], 2, poTrap.M)

# computation of the Jacobian at out_pof
_J1 = poTrap(Val(:JacFullSparse), outpo_f, par_hopf)
_Jfd = ForwardDiff.jacobian(z-> poTrap(z,par_hopf), outpo_f)

# test of the jacobian againt automatic differentiation
@test norm(_Jfd - Array(_J1), Inf) < 1e-7
####################################################################################################
# test PeriodicUtils
BK.amplitude(rand(10,10),3)
BK.amplitude(rand(101), 4, 25)
BK.maximumPOTrap(rand(101), 4, 25)
####################################################################################################
# tests for constructor of Floquet routines
BK.checkFloquetOptions(EigArpack())
BK.checkFloquetOptions(EigArnoldiMethod())
BK.checkFloquetOptions(EigKrylovKit())
FloquetQaD(EigKrylovKit()) |> FloquetQaD


text:
#! /usr/bin/env julia

module problem101

function polyval(coefficients :: Vector, x :: Number)
  powers = [x^p for p = 0:length(coefficients)-1]
  return sum(coefficients .* powers)
end

function polyfit(x :: Vector, y :: Vector)
  return polyfit(x, y, length(x))
end

function polyfit(x :: Vector, y :: Vector, n :: Int)
  A = [xi^p for xi = x, p = 0:(n-1)]
  return A \ y
end

function bad_opt_poly(poly, n)
  xs = [1:n]
  ys = [polyval(poly, x) for x = xs]
  return polyfit(xs, ys)
end

function bad_opt_polys(poly)
  return {bad_opt_poly(poly, n) for n = 1:length(poly)-1}
end

function first_incorrect_terms(polys)
  return [polyval(p, length(p)+1) for p = polys]
end

function sum_of_fits(poly)
  return sum(first_incorrect_terms(bad_opt_polys(poly)))
end

using Base.Test

@test polyval([1, 2, 3], 5) == 1 +  2*5 + 3*5*5
@test polyfit([1, 2, 3], [1, 8, 27]) == [6, -11, 6]
@test first_incorrect_terms(bad_opt_polys([0, 0, 0, 1])) == [1, 15, 58]
@test sum_of_fits([0, 0, 0, 1]) == 74

end

println(int(problem101.sum_of_fits([1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1])))

max_stars_repo_name: mbuhot/mbuhot-euler-solutions

max_issues_repo_name: mbuhot/mbuhot-euler-solutions

max_forks_repo_name: mbuhot/mbuhot-euler-solutions


text:
function _TestOut(model::Function, startp::AbstractVector, xlen::Int; max::Int=100)
    if MaximalNumberOfArguments(model) == 2
        model((xlen < 2 ? rand() : rand(xlen)), startp)
    else
        Res = fill(-Inf, max)
        model(Res, (xlen < 2 ? rand() : rand(xlen)), startp)
        Res[1:(findfirst(isinf, Res) - 1)]
    end
end

# Callback triggers when Boundaries is `true`.
"""
    ModelMap(Map::Function, InDomain::Union{Nothing,Function}, Domain::HyperCube)
    ModelMap(Map::Function, InDomain::Function, xyp::Tuple{Int,Int,Int})
A container which stores additional information about a model map, in particular its domain of validity.
`Map` is the actual map `(x,θ) -> model(x,θ)`. `Domain` is a `HyperCube` which allows one to roughly specify the ranges of the various parameters.
For more complicated boundary constraints, scalar function `InDomain` can be specified, which should be strictly positive on the valid parameter domain.
!!! note
    A `Bool`-valued function which returns `true` in the valid domain also fits this description, which allows one to easily combine multiple constraints.
    Providing this information about the domain can be advantageous in the optimization process for complicated models.
"""
struct ModelMap
    Map::Function
    InDomain::Union{Nothing,Function}
    Domain::Cuboid
    xyp::Tuple{Int,Int,Int}
    pnames::AbstractVector{<:String}
    StaticOutput::Val
    inplace::Val
    CustomEmbedding::Val
    # Given: Bool-valued domain function
    function ModelMap(model::Function, InDomain::Function, xyp::Tuple{Int,Int,Int}; pnames::Union{AbstractVector{<:String},Bool}=false)
        ModelMap(model, InDomain, nothing, xyp; pnames=pnames)
    end
    # Given: HyperCube
    function ModelMap(model::Function, Domain::Cuboid, xyp::Union{Tuple{Int,Int,Int},Bool}=false; pnames::Union{AbstractVector{<:String},Bool}=false)
        xyp isa Bool ? ModelMap(model, nothing, Domain; pnames=pnames) : ModelMap(model, nothing, Domain, xyp; pnames=pnames)
    end
    # Given: Function only (potentially) -> Find xyp
    function ModelMap(model::Function, InDomain::Union{Nothing,Function}=nothing, Domain::Union{Cuboid,Nothing}=nothing; pnames::Union{AbstractVector{<:String},Bool}=false)
        startp = isnothing(Domain) ? GetStartP(GetArgSize(model)[2]) : ElaborateGetStartP(Domain, InDomain)
        xlen = MaximalNumberOfArguments(model) > 2 ? GetArgLength((Res,x)->model(Res,x,startp)) : GetArgLength(x->model(x,startp))
        testout = _TestOut(model, startp, xlen)
        ModelMap(model, InDomain, Domain, (xlen, size(testout,1), length(startp)); pnames=pnames)
    end
    function ModelMap(model::Function, InDomain::Union{Nothing,Function}, Domain::Union{Cuboid,Nothing}, xyp::Tuple{Int,Int,Int}; pnames::Union{AbstractVector{<:String},Bool}=false)
        pnames = typeof(pnames) == Bool ? CreateSymbolNames(xyp[3],"θ") : pnames
        startp = isnothing(Domain) ? GetStartP(xyp[3]) : ElaborateGetStartP(Domain, InDomain)
        testout = _TestOut(model, startp, xyp[1])
        StaticOutput = testout isa SVector
        Inplace = MaximalNumberOfArguments(model) > 2
        ModelMap(model, InDomain, Domain, xyp, pnames, Val(StaticOutput), Val(Inplace), Val(false))
    end
    "Construct new ModelMap from function `F` with data from `M`."
    ModelMap(F::Function, M::ModelMap; inplace::Bool=isinplace(M)) = ModelMap(F, M.InDomain, M.Domain, M.xyp, M.pnames, M.StaticOutput, Val(inplace), M.CustomEmbedding)
    # Careful with inheriting CustomEmbedding to the Jacobian! For automatically generated dmodels (symbolic or autodiff) it should be OFF!
    function ModelMap(Map::Function, InDomain::Union{Nothing,Function}, Domain::Union{Cuboid,Nothing}, xyp::Tuple{Int,Int,Int},
                        pnames::AbstractVector{String}, StaticOutput::Val, inplace::Val=Val(false), CustomEmbedding::Val=Val(false))
        isnothing(Domain) && (Domain = FullDomain(xyp[3], 1e5))
        InDomain isa Function && (@assert InDomain(Center(Domain)) isa Number "InDomain function must yield a scalar value, got $(typeof(InDomain(Center(Domain)))) at $(Center(Domain)).")
        new(Map, InDomain, Domain, xyp, pnames, StaticOutput, inplace, CustomEmbedding)
    end
end
(M::ModelMap)(x, θ::AbstractVector{<:Number}; kwargs...) = M.Map(x, θ; kwargs...)
ModelOrFunction = Union{Function,ModelMap}

# For SciMLBase.remake
ModelMap(;
Map::Function=x->Inf,
InDomain::Union{Nothing,Function}=nothing,
Domain::Union{Cuboid,Nothing}=nothing,
xyp::Tuple{Int,Int,Int}=(1,1,1),
pnames::AbstractVector{String}=["θ"],
StaticOutput::Val=Val(true),
inplace::Val=Val(true),
CustomEmbedding::Val=Val(true)) = ModelMap(Map, InDomain, Domain, xyp, pnames, StaticOutput, inplace, CustomEmbedding)

function InformNames(M::ModelMap, pnames::AbstractVector{String})
    @assert length(pnames) == M.xyp[3]
    ModelMap(M.Map, M.InDomain, M.Domain, M.xyp, pnames, M.StaticOutput, M.inplace, M.CustomEmbedding)
end

pnames(M::ModelMap) = M.pnames
Domain(M::ModelMap) = M.Domain
isinplace(M::ModelMap) = ValToBool(M.inplace)
iscustom(M::ModelMap) = ValToBool(M.CustomEmbedding)

IsInDomain(M::ModelMap) = θ::AbstractVector -> IsInDomain(M, θ)
IsInDomain(M::ModelMap, θ::AbstractVector) = IsInDomain(M.InDomain, M.Domain, θ)
IsInDomain(InDomain::Union{Nothing,Function}, Domain::Union{Nothing,Cuboid}, θ::AbstractVector) = (_TestInDomain(InDomain, θ) && _TestDomain(Domain, θ))

# Eval InDomain function
_TestInDomain(::Nothing, θ::AbstractVector) = true
_TestInDomain(InDomain::Function, θ::AbstractVector) = InDomain(θ) > 0
# Eval Domain HyperCube
_TestDomain(::Nothing, θ::AbstractVector) = true       # Excluded
_TestDomain(Domain::Cuboid, θ::AbstractVector) = θ ∈ Domain

MakeCustom(F::Function, Domain::Union{Cuboid,Bool,Nothing}=nothing) = Domain isa Cuboid ? MakeCustom(ModelMap(F, Domain)) : MakeCustom(ModelMap(F))
function MakeCustom(M::ModelMap)
    if iscustom(M)
        println("Map already uses custom embedding.")
        return M
    else
        return ModelMap(M.Map, M.InDomain, M.Domain, M.xyp, M.pnames, M.StaticOutput, M.inplace, Val(true))
    end
end
function MakeNonCustom(M::ModelMap)
    if !iscustom(M)
        println("Map already not using custom embedding.")
        return M
    else
        return ModelMap(M.Map, M.InDomain, M.Domain, M.xyp, M.pnames, M.StaticOutput, M.inplace, Val(false))
    end
end

function ModelMap(F::Nothing, M::ModelMap)
    println("ModelMap: Got nothing instead of function to build new ModelMap")
    nothing
end
function CreateSymbolNames(n::Int, base::String="θ")
    n == 1 && return [base]
    D = Dict(string.(0:9) .=> ["₀","₁","₂","₃","₄","₅","₆","₇","₈","₉"])
    base .* [prod(get(D,"$x","Q") for x in string(digit)) for digit in 1:n]
end

pdim(DS::AbstractDataSet, model::ModelMap)::Int = model.xyp[3]

function ModelMappize(DM::AbstractDataModel)
    NewMod = Predictor(DM) isa ModelMap ? Predictor(DM) : ModelMap(Predictor(DM))
    NewdMod = dPredictor(DM) isa ModelMap ? dPredictor(DM) : ModelMap(dPredictor(DM))
    DataModel(Data(DM), NewMod, NewdMod, MLE(DM), LogLikeMLE(DM))
end

"""
Only works for `DataSet` and `DataSetExact` but will output wrong order of components for `CompositeDataSet`!
"""
function ConcatenateModels(Mods::AbstractVector{<:ModelMap})
    @assert ConsistentElDims((x->x.xyp[1]).(Mods)) > 0 && ConsistentElDims((x->x.xyp[3]).(Mods)) > 0
    if Mods[1].xyp[1] == 1
        function ConcatenatedModel(x::Number, θ::AbstractVector{<:Number}; kwargs...)
            map(model->model(x, θ; kwargs...), Mods) |> Reduction
        end
        EbdMap(model::Function, θ::AbstractVector, woundX::AbstractVector,custom::Val{false}; kwargs...) = Reduction(map(x->model(x,θ; kwargs...), woundX))
        EbdMap(model::Function, θ::AbstractVector, woundX::AbstractVector,custom::Val{true}; kwargs...) = model(woundX, θ; kwargs...)
        function ConcatenatedModel(X::AbstractVector{<:Number}, θ::AbstractVector{<:Number}; kwargs...)
            if any(iscustom, Mods)
                Res = if any(m->m.xyp[2]>1, Mods)
                    map(m->Windup(EbdMap(m.Map, θ, X, m.CustomEmbedding; kwargs...), m.xyp[2]), Mods)
                    # map(m->Windup(EmbeddingMap(DS, m, θ, X), m.xyp[2]), Mods)
                else
                    map(m->EbdMap(m.Map, θ, X, m.CustomEmbedding; kwargs...), Mods)
                    # map(m->EmbeddingMap(DS, m, θ, X), Mods)
                end
                return zip(Res...) |> Iterators.flatten |> collect |> Reduction
            else
                return map(z->ConcatenatedModel(z, θ; kwargs...), X) |> Reduction
            end
        end
        return ModelMap(ConcatenatedModel, reduce(union, (z->z.Domain).(Mods)), (Mods[1].xyp[1], sum((q->q.xyp[2]).(Mods)), Mods[1].xyp[3])) |> MakeCustom
    else
        function NConcatenatedModel(x::AbstractVector{<:Number}, θ::AbstractVector{<:Number}; kwargs...)
            map(model->model(x, θ; kwargs...), Mods) |> Reduction
        end
        function NConcatenatedModel(X::AbstractVector{<:AbstractVector{<:Number}}, θ::AbstractVector{<:Number}; kwargs...)
            if any(iscustom, Mods)
                Res = if any(m->m.xyp[2]>1, Mods)
                    map(m->Windup(EmbeddingMap(DS, m, θ, X), m.xyp[2]), Mods)
                else
                    map(m->EmbeddingMap(DS, m, θ, X), Mods)
                end
                return zip(Res...) |> Iterators.flatten |> collect |> Reduction
            else
                return map(z->NConcatenatedModel(z, θ; kwargs...), X) |> Reduction
            end
        end
        return ModelMap(NConcatenatedModel, reduce(union, (z->z.Domain).(Mods)), (Mods[1].xyp[1], sum((q->q.xyp[2]).(Mods)), Mods[1].xyp[3])) |> MakeCustom
    end
end

_Apply(x::AbstractVector{<:Number}, Componentwise::Function, idxs::BoolVector) = [(idxs[i] ? Componentwise(x[i]) : x[i]) for i in eachindex(idxs)]
_ApplyFull(x::AbstractVector{<:Number}, Vectorial::Function) = Vectorial(x)

MonotoneIncreasing(F::Function, Interval::Tuple{Number,Number})::Bool = Monotonicity(F, Interval) == :increasing
MonotoneDecreasing(F::Function, Interval::Tuple{Number,Number})::Bool = Monotonicity(F, Interval) == :decreasing
function Monotonicity(F::Function, Interval::Tuple{Number,Number})
    derivs = map(x->ForwardDiff.derivative(F, x), range(Interval[1], Interval[2]; length=200))
    all(x-> x≥0., derivs) && return :increasing
    all(x-> x≤0., derivs) && return :decreasing
    :neither
end

Transform(model::Function, idxs::BoolVector, Transform::Function, InverseTransform::Function=x->invert(Transform,x)) = _Transform(model, idxs, Transform, InverseTransform)

# Try to do a bit of inference for the new domain here!
function Transform(M::ModelMap, idxs::BoolVector, Transform::Function, InverseTransform::Function=x->invert(Transform,x))
    TransformedDomain = M.InDomain isa Function ? (θ::AbstractVector{<:Number} -> M.InDomain(_Apply(θ, Transform, idxs))) : nothing
    mono = Monotonicity(Transform, (1e-12,50.))
    NewCube = if mono == :increasing
        HyperCube(_Apply(M.Domain.L, InverseTransform, idxs), _Apply(M.Domain.U, InverseTransform, idxs))
    elseif mono == :decreasing
        println("Detected monotone decreasing transformation.")
        HyperCube(_Apply(M.Domain.U, InverseTransform, idxs), _Apply(M.Domain.L, InverseTransform, idxs))
    else
        @warn "Transformation does not appear to be monotone. Unable to infer new Domain."
        FullDomain(length(idxs))
    end
    ModelMap(_Transform(M.Map, idxs, Transform, InverseTransform), TransformedDomain, NewCube,
                        M.xyp, M.pnames, M.StaticOutput, M.inplace, M.CustomEmbedding)
end
# function Transform(M::ModelMap, Transform::Function, InverseTransform::Function=x->invert(Transform,x))
#     Transform(M, trues(M.xyp[3]), Transform, InverseTransform)
# end

function _Transform(F::Function, idxs::BoolVector, Transform::Function, InverseTransform::Function)
    function TransformedModel(x::Union{Number, AbstractVector{<:Number}}, θ::AbstractVector{<:Number}; kwargs...)
        F(x, _Apply(θ, Transform, idxs); kwargs...)
    end
end

"""
    Transform(DM::AbstractDataModel, F::Function, idxs=trues(pdim(DM))) -> DataModel
    Transform(model::Function, idxs, F::Function) -> Function
Transforms the parameters of the model by the given scalar function `F` such that `newmodel(x, θ) = oldmodel(x, F.(θ))`.
By providing `idxs`, one may restrict the application of the function `F` to specific parameter components.
"""
function Transform(DM::AbstractDataModel, F::Function, idxs::BoolVector=trues(pdim(DM)); kwargs...)
    @assert length(idxs) == pdim(DM)
    sum(idxs) == 0 && return DM
    DataModel(Data(DM), Transform(Predictor(DM), idxs, F), _Apply(MLE(DM), x->invert(F,x), idxs); kwargs...)
end
function Transform(DM::AbstractDataModel, F::Function, inverseF::Function, idxs::BoolVector=trues(pdim(DM)); kwargs...)
    @assert length(idxs) == pdim(DM)
    sum(idxs) == 0 && return DM
    DataModel(Data(DM), Transform(Predictor(DM), idxs, F, inverseF), _Apply(MLE(DM), inverseF, idxs); kwargs...)
end

LogTransform(M::ModelOrFunction, idxs::BoolVector=(M isa ModelMap ? trues(M.xyp[3]) : trues(GetArgSize(M)[2]))) = Transform(M, idxs, log, exp)
LogTransform(DM::AbstractDataModel, idxs::BoolVector=trues(pdim(DM)); kwargs...) = Transform(DM, log, exp, idxs; kwargs...)

ExpTransform(M::ModelOrFunction, idxs::BoolVector=(M isa ModelMap ? trues(M.xyp[3]) : trues(GetArgSize(M)[2]))) = Transform(M, idxs, exp, log)
ExpTransform(DM::AbstractDataModel, idxs::BoolVector=trues(pdim(DM)); kwargs...) = Transform(DM, exp, log, idxs; kwargs...)

Log10Transform(M::ModelOrFunction, idxs::BoolVector=(M isa ModelMap ? trues(M.xyp[3]) : trues(GetArgSize(M)[2]))) = Transform(M, idxs, log10, x->10^x)
Log10Transform(DM::AbstractDataModel, idxs::BoolVector=trues(pdim(DM)); kwargs...) = Transform(DM, log10, x->10^x, idxs; kwargs...)

Power10Transform(M::ModelOrFunction, idxs::BoolVector=(M isa ModelMap ? trues(M.xyp[3]) : trues(GetArgSize(M)[2]))) = Transform(M, idxs, x->10^x, log10)
Power10Transform(DM::AbstractDataModel, idxs::BoolVector=trues(pdim(DM)); kwargs...) = Transform(DM, x->10^x, log10, idxs; kwargs...)

ReflectionTransform(M::ModelOrFunction, idxs::BoolVector=(M isa ModelMap ? trues(M.xyp[3]) : trues(GetArgSize(M)[2]))) = Transform(M, idxs, x-> -x, x-> -x)
ReflectionTransform(DM::AbstractDataModel, idxs::BoolVector=trues(pdim(DM)); kwargs...) = Transform(DM, x-> -x, x-> -x, idxs; kwargs...)

ScaleTransform(M::ModelOrFunction, factor::Number, idxs::BoolVector=(M isa ModelMap ? trues(M.xyp[3]) : trues(GetArgSize(M)[2]))) = Transform(M, idxs, x->factor*x, x->x/factor)
ScaleTransform(DM::AbstractDataModel, factor::Number, idxs::BoolVector=trues(pdim(DM)); kwargs...) = Transform(DM, x->factor*x, x->x/factor, idxs; kwargs...)

function TranslationTransform(F::Function, v::AbstractVector{<:Number})
    TranslatedModel(x, θ::AbstractVector{<:Number}; kwargs...) = F(x, θ + v; kwargs...)
end
function TranslationTransform(M::ModelMap, v::AbstractVector{<:Number})
    @assert length(M.Domain) == length(v)
    ModelMap(TranslationTransform(M.Map, v), (M.InDomain isa Function ? (θ->M.InDomain(θ + v)) : nothing), TranslateCube(M.Domain, -v), M.xyp, M.pnames, M.StaticOutput,
                                    M.inplace, M.CustomEmbedding)
end
function TranslationTransform(DM::AbstractDataModel, v::AbstractVector{<:Number}; kwargs...)
    @assert pdim(DM) == length(v)
    DataModel(Data(DM), TranslationTransform(Predictor(DM), v), MLE(DM)-v; kwargs...)
end

function LinearTransform(F::Function, A::AbstractMatrix{<:Number})
    TransformedModel(x, θ::AbstractVector{<:Number}; kwargs...) = F(x, A*θ; kwargs...)
end
function LinearTransform(M::ModelMap, A::AbstractMatrix{<:Number})
    @assert length(M.Domain) == size(A,1) == size(A,2)
    Ainv = inv(A)
    ModelMap(LinearTransform(M.Map, A), (M.InDomain isa Function ? (θ->M.InDomain(A*θ)) : nothing), HyperCube(Ainv * M.Domain.L, Ainv * M.Domain.U),
                    M.xyp, M.pnames, M.StaticOutput, M.inplace, M.CustomEmbedding)
end
function LinearTransform(DM::AbstractDataModel, A::AbstractMatrix{<:Number}; kwargs...)
    @assert pdim(DM) == size(A,1) == size(A,2)
    DataModel(Data(DM), LinearTransform(Predictor(DM), A), inv(A)*MLE(DM); kwargs...)
end

function AffineTransform(F::Function, A::AbstractMatrix{<:Number}, v::AbstractVector{<:Number})
    @assert size(A,1) == size(A,2) == length(v)
    TranslatedModel(x, θ::AbstractVector{<:Number}; kwargs...) = F(x, A*θ + v; kwargs...)
end
function AffineTransform(M::ModelMap, A::AbstractMatrix{<:Number}, v::AbstractVector{<:Number})
    @assert length(M.Domain) == size(A,1) == size(A,2) == length(v)
    Ainv = inv(A)
    ModelMap(AffineTransform(M.Map, A, v), (M.InDomain isa Function ? (θ->M.InDomain(A*θ+v)) : nothing), HyperCube(Ainv*(M.Domain.L-v), Ainv*(M.Domain.U-v)),
                    M.xyp, M.pnames, M.StaticOutput, M.inplace, M.CustomEmbedding)
end
function AffineTransform(DM::AbstractDataModel, A::AbstractMatrix{<:Number}, v::AbstractVector{<:Number}; kwargs...)
    @assert pdim(DM) == size(A,1) == size(A,2) == length(v)
    Ainv = inv(A)
    DataModel(Data(DM), AffineTransform(Predictor(DM), A, v), Ainv*(MLE(DM)-v); kwargs...)
end

_GetDecorrelationTransform(DM::AbstractDataModel) = (_GetDecorrelationTransform(FisherMetric(DM, MLE(DM))), MLE(DM))
_GetDecorrelationTransform(M::AbstractMatrix) = cholesky(Symmetric(inv(M))).L
LinearDecorrelation(DM::AbstractDataModel; kwargs...) = ((M,X) =_GetDecorrelationTransform(DM); AffineTransform(DM, M, X; kwargs...))

function DecorrelationTransforms(DM::AbstractDataModel)
    M, X = _GetDecorrelationTransform(DM);      iM = inv(M)
    ForwardTransform(x::AbstractVector) = M * x + X
    InvTransform(x::AbstractVector) = iM * (x - X)
    ForwardTransform, InvTransform
end

"""
    EmbedModelVia(model, F::Function; Domain::HyperCube=FullDomain(GetArgLength(F))) -> Union{Function,ModelMap}
Transforms a model function via `newmodel(x, θ) = oldmodel(x, F(θ))`.
A `Domain` for the new model can optionally be specified for `ModelMap`s.
"""
function EmbedModelVia(model::Function, F::Function; Kwargs...)
    EmbeddedModel(x, θ; kwargs...) = model(x, F(θ); kwargs...)
end
function EmbedModelVia(M::ModelMap, F::Function; Domain::Union{Nothing,HyperCube}=nothing)
    if isnothing(Domain)
        Domain = FullDomain(GetArgLength(F))
        @warn "Cannot infer new Domain HyperCube for general embeddings, using $Domain."
    end
    ModelMap(EmbedModelVia(M.Map, F), (M.InDomain isa Function ? (M.InDomain∘F) : nothing), Domain, (M.xyp[1], M.xyp[2], length(Domain)), CreateSymbolNames(length(Domain), "θ"), M.StaticOutput, M.inplace, M.CustomEmbedding)
end
function EmbedDModelVia(dmodel::Function, F::Function; ADmode::Union{Symbol,Val}=Val(:ForwardDiff), Kwargs...)
    Jac = GetJac(ADmode,F)
    EmbeddedJacobian(x, θ; kwargs...) = dmodel(x, F(θ); kwargs...) * Jac(θ)
end
function EmbedDModelVia(dM::ModelMap, F::Function; Domain::HyperCube=FullDomain(GetArgLength(F)))
    ModelMap(EmbedDModelVia(dM.Map, F), (dM.InDomain isa Function ? (dM.InDomain∘F) : nothing), Domain, (dM.xyp[1], dM.xyp[2], length(Domain)), CreateSymbolNames(length(Domain), "θ"), dM.StaticOutput, dM.inplace, dM.CustomEmbedding)
end

"""
    Embedding(DM::AbstractDataModel, F::Function, start::AbstractVector; Domain::HyperCube=FullDomain(length(start))) -> DataModel
Transforms a model function via `newmodel(x, θ) = oldmodel(x, F(θ))` and returns the associated `DataModel`.
An initial parameter configuration `start` as well as a `Domain` can optionally be passed to the `DataModel` constructor.
"""
function Embedding(DM::AbstractDataModel, F::Function, start::AbstractVector{<:Number}=GetStartP(GetArgLength(F)); Domain::HyperCube=FullDomain(length(start)), kwargs...)
    DataModel(Data(DM), EmbedModelVia(Predictor(DM), F; Domain=Domain), EmbedDModelVia(dPredictor(DM), F; Domain=Domain), start; kwargs...)
end

ExpTransform(Sys::ODESystem, idxs::AbstractVector{<:Bool}=trues(length(parameters(Sys))); kwargs...) = SystemTransform(Sys, exp, idxs; kwargs...)
LogTransform(Sys::ODESystem, idxs::AbstractVector{<:Bool}=trues(length(parameters(Sys))); kwargs...) = SystemTransform(Sys, log, idxs; kwargs...)
Power10Transform(Sys::ODESystem, idxs::AbstractVector{<:Bool}=trues(length(parameters(Sys))); kwargs...) = SystemTransform(Sys, x->(10.0)^x, idxs; kwargs...)
Log10Transform(Sys::ODESystem, idxs::AbstractVector{<:Bool}=trues(length(parameters(Sys))); kwargs...) = SystemTransform(Sys, log10, idxs; kwargs...)

"""
    SystemTransform(Sys::ODESystem, F::Function, idxs::AbstractVector{<:Bool}=trues(length(parameters(Sys)))) -> ODESystem
Transforms the parameters of a `ODESystem` according to `F`.
"""
function SystemTransform(Sys::ODESystem, F::Function, idxs::AbstractVector{<:Bool}=trues(length(parameters(Sys))))
    SubstDict = Dict(parameters(Sys) .=> [(idxs[i] ? F(x) : x) for (i,x) in enumerate(parameters(Sys))])
    NewEqs = [(equations(Sys)[i].lhs ~ substitute(equations(Sys)[i].rhs, SubstDict)) for i in 1:length(equations(Sys))]
    ODESystem(NewEqs, independent_variables(Sys)[1], states(Sys), parameters(Sys); name=nameof(Sys))
end

LinearModel(x::Union{Number,AbstractVector{<:Number}}, θ::AbstractVector{<:Number}) = dot(θ[1:end-1], x) + θ[end]
QuadraticModel(x::Union{Number,AbstractVector{<:Number}}, θ::AbstractVector{<:Number}) = dot(θ[1:Int((end-1)/2)], x.^2) + dot(θ[Int((end-1)/2)+1:end-1], x) + θ[end]
ExponentialModel(x::Union{Number,AbstractVector{<:Number}}, θ::AbstractVector{<:Number}) = exp(LinearModel(x,θ))
SumExponentialsModel(x::Union{Number,AbstractVector{<:Number}}, θ::AbstractVector{<:Number}) = sum(exp.(θ .* x))

function PolynomialModel(degree::Int)
    Polynomial(x::Number, θ::AbstractVector{<:Number}) = sum(θ[i] * x^(i-1) for i in 1:(degree+1))
end


text:
struct AffineMap

    Mat :: SparseMatrixCSC{Float64, Int64}
    f :: Vector{Float64}

end

function Base.:*(A::AffineMap, v::Vector)
    A.Mat*v + A.f
end

function Base.:*(M::SparseMatrixCSC{Float64, Int64}, A::AffineMap)
    AffineMap(M*A.Mat, M*A.f)
end

function Base.:*(A::AffineMap, M::SparseMatrixCSC{Float64, Int64})
    AffineMap(A.Mat*M, A.f)
end

function Base.:*(A1::AffineMap, A2::AffineMap)
    Mat = A1.Mat*A2.Mat
    f = A1.Mat*A2.f + A1.f
    AffineMap(Mat, f)
end

mutable struct LinearProblem

    A :: SparseMatrixCSC{Float64, Int64}
    Adiag :: Vector{Float64}
    b :: Vector{Float64}
    x :: Vector{Float64}
    r :: Vector{Float64}

end
function LinearProblem(A::SparseMatrixCSC, x::Vector{Float64}, b::Vector{Float64})
    Adiag = diag(A)
    r = b - A*x
    return LinearProblem(A, Adiag, b, x, r)
end

function compute_residual!(L::LinearProblem)
    L.r[:] = L.b - L.A*L.x
end

function coarsen_problem!(L::LinearProblem, Restr::AffineMap, Interp::AffineMap)

    lhs = Restr * (L.A*Interp)
    rhs = Restr * L.b

    A = lhs.Mat
    b = rhs - lhs.f
    x = zeros(length(b))

    return LinearProblem(A, x, b)
end

function jsmooth!(L::LinearProblem, n::Int64)

    for _=1:n
        L.x[:] += L.r ./ L.Adiag
        compute_residual!(L)
    end
end


text:
using Dates
using Random: randstring

"""
    unshortlink(url)

return unshorten url or the url if it is not a short link
"""
function unshortlink(url; kw...)
    rq = HTTP.request("HEAD", url; redirect=false, status_exception=false, kw...)
    while rq.status ÷ 100 == 3
        url = HTTP.header(rq, "Location")
        rq = HTTP.request("HEAD", url; redirect=false, status_exception=false, kw...)
    end
    url
end

isgooglesheet(url) = occursin("docs.google.com/spreadsheets", url)
isgoogledrive(url) = occursin("docs.google.com", url)

function googlesheet_handler(url; format=:csv)
    link, expo = splitdir(url)
    if startswith(expo, "edit") || expo == ""
        url = link * "/export?format=$format"
    elseif startswith(expo, "export")
        url = replace(url, r"format=([a-zA-Z]*)(.*)"=>SubstitutionString("format=$format\\2"))
    end
    url
end

function maybegoogle_download(url, localdir)
    long_url = unshortlink(url)
    if isgooglesheet(long_url)
        long_url = googlesheet_handler(long_url)
    end

    if isgoogledrive(long_url)
        download_gdrive(long_url, localdir)
    else
        DataDeps.fetch_http(long_url, localdir)
    end
end

function find_gcode(ckj)
    for cookie ∈ ckj
        if match(r"download_warning_", cookie.name) !== nothing
            return cookie.value
        end
    end

    nothing
end

function download_gdrive(url, localdir; retry=true, retries=4)
    gcode = nothing
    try_time = 0
    ckjar = Dict{String, Set{HTTP.Cookies.Cookie}}()
    while isnothing(gcode) && retry && try_time < retries
        if try_time > 0
            @info "retrying..."
        end
        try_time += 1
        rq = HTTP.request("HEAD", url; cookies=true, cookiejar=ckjar)
        ckj = ckjar["docs.google.com"]
        gcode = find_gcode(ckj)
        if isnothing(gcode)
            @warn "gcode not found." rq
            sleep(3)
        end
    end

    if isnothing(gcode)
        error("download failed")
    end

    format_progress(x) = round(x, digits=4)
    format_bytes(x) = !isfinite(x) ? "∞ B" : Base.format_bytes(x)
    format_seconds(x) = "$(round(x; digits=2)) s"
    format_bytes_per_second(x) = format_bytes(x) * "/s"

    local filepath
    # newurl = unshortlink("$url&confirm=$gcode"; cookies=true, cookiejar=ckjar)
    newurl = "$url&confirm=$gcode"

    #part of codes are from https://github.com/JuliaWeb/HTTP.jl/blob/master/src/download.jl
    HTTP.open("GET", newurl, ["Range"=>"bytes=0-"]; cookies=true, cookiejar=ckjar, redirect_limit=10) do stream
        resp = HTTP.startread(stream)
        hcd = HTTP.header(resp, "Content-Disposition")
        isempty(hcd) && return

        m = match(r"filename=\\\"(.*)\\\"", hcd)
        if m === nothing
            filename = "gdrive_downloaded-$(randstring())"
        else
            filename = m.captures[]
        end

        filepath = joinpath(localdir, filename)

        total_bytes = tryparse(Float64, split(HTTP.header(resp, "Content-Range"), '/')[end])
        total_bytes === nothing && (total_bytes = NaN)
        downloaded_bytes = 0
        start_time = now()
        prev_time = now()
        period = DataDeps.progress_update_period()

        function report_callback()
            prev_time = now()
            taken_time = (prev_time - start_time).value / 1000 # in seconds
            average_speed = downloaded_bytes / taken_time
            remaining_bytes = total_bytes - downloaded_bytes
            remaining_time = remaining_bytes / average_speed
            completion_progress = downloaded_bytes / total_bytes

            @info("Downloading",
                  source=url,
                  dest = filepath,
                  progress = completion_progress |> format_progress,
                  time_taken = taken_time |> format_seconds,
                  time_remaining = remaining_time |> format_seconds,
                  average_speed = average_speed |> format_bytes_per_second,
                  downloaded = downloaded_bytes |> format_bytes,
                  remaining = remaining_bytes |> format_bytes,
                  total = total_bytes |> format_bytes,
                  )
        end

        Base.open(filepath, "w") do fh
            while(!eof(stream))
                downloaded_bytes += write(fh, readavailable(stream))
                if !isinf(period)
                  if now() - prev_time > Millisecond(1000*period)
                    report_callback()
                  end
                end
            end
        end
        report_callback()
    end
    filepath
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMPlates.jl/blob/master/LICENSE

struct MindlinPlate <: FieldProblem end

function FEMBase.get_unknown_field_name(::Problem{MindlinPlate})
    return "displacement"
end

function FEMBase.assemble_elements!(problem::Problem{MindlinPlate},
                                    assembly::Assembly,
                                    elements::Vector{Element{Quad4}},
                                    time::Float64)

    for element in elements
        ndofs = 4*3
        Ke = zeros(ndofs, ndofs)
        fe = zeros(ndofs)
        for ip in get_integration_points(element)
            detJ = element(ip, time, Val{:detJ})
            E = element("youngs modulus", ip, time)
            t = element("thickness", ip, time)
            nu = element("poissons ratio", ip, time)
            gamma = 5/6
            if haskey(element, "shear correction factor")
                gamma = element("shear correction factor", ip, time)
            end
            # http://what-when-how.com/the-finite-element-method/fem-for-plates-and-shells-finite-element-method-part-1/
            G = E/(2.0 * (1.0 + nu))
            Db = E*t^3/(12*(1-nu^2)) * [
                1.0  nu  0.0
                nu  1.0  0.0
                0.0 0.0 (1-nu)/2]
            Ds = G*t*gamma*[1.0 0.0; 0.0 1.0]
            B_I = zeros(3, 12)
            N = element(ip, time)
            Grad = element(ip, time, Val{:Grad})
            for j=1:4
                B_I[1, (j-1)*3+2] = Grad[1,j]
                B_I[2, (j-1)*3+3] = Grad[2,j]
                B_I[3, (j-1)*3+2] = Grad[2,j]
                B_I[3, (j-1)*3+3] = Grad[1,j]
            end
            B_O = zeros(2, 12)
            for j=1:4
                B_O[1, (j-1)*3+1] = Grad[1,j]
                B_O[1, (j-1)*3+2] = N[j]
                B_O[2, (j-1)*3+1] = Grad[2,j]
                B_O[2, (j-1)*3+3] = N[j]
            end
            z = zeros(3, 2)
            D = [Db z; z' Ds]
            K1 = ip.weight * B_I'*Db*B_I * detJ
            K2 = ip.weight * B_O'*Ds*B_O * detJ
            Ke += (K1 + K2)
            if haskey(element, "distributed load")
                p = element("distributed load", ip, time)
                fe[1:3:end] += ip.weight * p * N' * detJ
            end
        end
        gdofs = get_gdofs(problem, element)
        add!(assembly.K, gdofs, gdofs, Ke)
        add!(assembly.f, gdofs, fe)
    end

end


text:
export delaunayWall;
#===============================================================================
#
#	src/deWall.jl
#
#	This File Contains:
#
#	 - delaunayWall(
#			P::Lar.Points,
#			ax = 1,
#			Pblack = Float64[],
#			AFL = Array{Int64,1}[],
#			tetraDict = DataStructures.Dict{Array{Int64,1},Array{Float64,1}}();
#			DEBUG = false
#		)::Lar.Cells
#
#	 - firstDeWallSimplex(
#			P::Lar.Points,
#			ax::Int64,
#			off::Float64;
#			DEBUG = false
#		)::Array{Int64,1}
#
#	 - findWallSimplex(
#			P::Lar.Points,
#			blackidx::Int64,
#			face::Array{Int64,1},
#			oppoint::Array{Float64,1};
#			DEBUG = false
#		)::Union{Array{Int64,1}, Nothing}
#
#	 - recursiveDelaunayWall(
# 			P::Lar.Points,
#			Pblack::Lar.Points,
# 			tetraDict::DataStructures.Dict{Array{Int64,1},Array{Float64,1}},
# 			AFL::Array{Array{Int64,1},1},
# 			ax::Int64,
# 			off::Float64,
# 			positive::Bool;
# 			DEBUG = false
# 		)::Lar.Cells
#
#	 - updateAFL!(
#			P::Lar.Points
#			new::Array{Int64,1}[],
#			AFLα = Array{Int64,1}[],
#			AFLplus = Array{Int64,1}[],
#			AFLminus = Array{Int64,1}[],
#			ax::Int64,
#			off::Float64;
# 			DEBUG = false
#		)::Bool
#
#	 - updatelist!(
#			list,
#			element
#		)::Bool
#
#	 - updateTetraDict!(
#			P::Lar.Points,
#			tetraDict::DataStructures.Dict{Array{Int64,1},Array{Float64,1}},
#			AFL::Array{Array{Int64,1},1},
#			σ::Array{Int64,1}
#		)::Nothing
#
===============================================================================#

"""
	delaunayWall(
		P::Lar.Points, ax = 1, Pblack::Float64[], AFL = Array{Int64,1}[],
		tetraDict = DataStructures.Dict{Array{Int64,1},Array{Float64,1}}();
		DEBUG = false
	)::Lar.Cells

Return the Delaunay Triangulation of sites `P` via Delaunay Wall algorithm.
The optional argument `ax` specify on wich axis it will build the Wall.
The optional argument `AFL` is used in recursive call.
If the keyword argument `DEBUG` is set to true than all the procedure is shown.

# Examples
```jldoctest

julia> P = [
 0.0 1.0 0.0 1.0 0.0 1.0 0.0 1.0
 0.0 0.0 1.0 1.0 0.0 0.0 1.0 1.0
 0.0 0.0 0.0 0.0 1.0 1.0 1.0 1.0
];

julia> DT = AlphaShapes.delaunayWall(P)
6-element Array{Array{Int64,1},1}:
 [1, 2, 4, 5]
 [1, 3, 4, 5]
 [2, 4, 5, 6]
 [3, 4, 5, 7]
 [4, 5, 6, 7]
 [4, 6, 7, 8]

```
"""

@timeit to "computeFirstSimplex" function computeFirstSimplex(P::Lar.Points,
	ax::Int64,
	off::Float64,
	tetraDict::Dict{Array{Int64,1},Array{Float64,1}},
	DT::Array{Array{Int64,1},1},
	DEBUG::Bool)::Array{Array{Int64,1},1}

	σ = sort(AlphaShapes.firstDeWallSimplex(P, ax, off, DEBUG = DEBUG))
	push!(DT, σ)
	AFL = @spawn AlphaShapes.simplexFaces(σ)
	AFL = fetch(AFL)

	AlphaShapes.updateTetraDict!(P, tetraDict, AFL, σ)
	return AFL
end

@timeit to "delaunayWall" function delaunayWall(
		P::Lar.Points,
		ax = 1,
		Pblack = Float64[],
		AFL = Array{Int64,1}[],
		tetraDict = DataStructures.Dict{Array{Int64,1},Array{Float64,1}}();
		DEBUG = false
	)::Lar.Cells

	if DEBUG @show "Delaunay Wall with parameters" P ax AFL tetraDict end

	# 0 - Data Reading and Container definition
	DT = Array{Int64,1}[]		# Delaunay Triangulation
	AFLα = Array{Int64,1}[]		# (d-1)faces intersecting the Wall
	AFLplus = Array{Int64,1}[]  # (d-1)faces in positive Wall half-space
	AFLminus = Array{Int64,1}[] # (d-1)faces in positive Wall half-space
	off =  AlphaShapes.findMedian(P, ax)

	if !isempty(Pblack) Pext = [P Pblack] else Pext = copy(P) end

	# 1 - Determine first simplex (if necessary)
	if isempty(AFL)
		@assert isempty(Pblack) "delaunayWall: If AFL is empty => Pblack must be"
		@assert isempty(tetraDict) "delaunayWall: If AFL is empty => tetraDict must be"
		AFL =  computeFirstSimplex(P,ax,off,tetraDict,DT,DEBUG)
	else
		@assert !isempty(Pblack) "delaunayWall: Data missing - Pblack"
		@assert !isempty(AFL) "delaunayWall: Data missing - AFL"
		@assert !isempty(tetraDict) "delaunayWall: Data missing - tetraDict"
	end

	# 2 - Build `AFL*` according to the axis `ax` with constant term `off`
	AlphaShapes.updateAFL!(
		P, AFL, AFLα, AFLplus, AFLminus, ax, off, DEBUG = DEBUG
	)

	# 4 - Build simplex Wall
	while !isempty(AFLα)
		# if face ∈ keys(tetraDict) oppoint = tetraDict[face]
		# else Pselection = setdiff([i for i = 1 : n], face) end

		σ =  AlphaShapes.findWallSimplex(
				Pext, AFLα[1], tetraDict[AFLα[1]], size(P, 2), DEBUG = DEBUG
			)
		if σ != nothing && σ ∉ DT
			push!(DT, σ)
			AFL = @spawn AlphaShapes.simplexFaces(σ)
			AFL = fetch(AFL)

			AlphaShapes.updateTetraDict!(P, tetraDict, AFL, σ)
			# Split σ's Faces according in semi-spaces
			AlphaShapes.updateAFL!(
				P, AFL, AFLα, AFLplus, AFLminus, ax, off, DEBUG=DEBUG
			)
		else
			@assert AlphaShapes.updatelist!(AFLα, AFLα[1]) == false "delaunayWall:
				Something unespected happends while removing a face."
		end
	end

	# 5 - Change the axis `ax` and repeat until there are no faces but exposed.
	#      A.K.A. Divide & Conquer phase.

	if !isempty(AFLminus)
		res = @spawn recursiveDelaunayWall(P, Pblack, tetraDict, AFLminus, ax, off, false; DEBUG = DEBUG)
		union!(DT, fetch(res))
	end

	if !isempty(AFLplus)
		res = @spawn recursiveDelaunayWall(P, Pblack, tetraDict, AFLplus, ax, off, true; DEBUG = DEBUG)
		union!(DT, fetch(res))
	end

	return DT
end

#-------------------------------------------------------------------------------
"""
	findWallSimplex(
		P::Lar.Points,
		face::Array{Int64,1}, oppoint::Array{Float64,1},
		blackidx = size(P, 2);
		DEBUG = false
	)::Union{Array{Int64,1}, Nothing}

Returns the simplex 'σ' build with `face` and a point from `P[:, 1:blackidx]`
such that it is in the opposite half plane of `oppoint`.
If such a simplex do not exists it returns `nothing`.
If `blackidx` is not specified all the points `P` are treaten as valid.
If the keyword argument `DEBUG` is set to true than all the procedure is shown.

# Examples
```jldoctest

julia> P = [
 0. 1. 0. 0. 2.;
 0. 0. 1. 0. 2.;
 0. 0. 0. 1. 2.
];

julia> newtetra = AlphaShapes.findWallSimplex(P,[2,3,4],[0., 0., 0.])
4-element Array{Int64,1}:
 2
 3
 4
 5

```
"""
@timeit to "findWallSimplex" function findWallSimplex(
		P::Lar.Points,
		face::Array{Int64,1},
		oppoint::Array{Float64,1},
		blackidx = size(P, 2);
		DEBUG = false
	)::Union{Array{Int64,1}, Nothing}

	if DEBUG @show "find Wall Simplex of" face oppoint end
	# Find the points in the halfspace defined by `face` that do not
	#  containsother the other point of the simplex.
	Pselection =
		@spawn AlphaShapes.oppositeHalfSpacePoints(P, P[:, face], oppoint)
	Pselection = fetch(Pselection)

	if DEBUG @show Pselection end

	# If there are no such points than the face is part of the convex hull.
	if isempty(Pselection)
		return nothing
	end

	# Find the Closest Point in the other halfspace with respect to σ
	#  according to dd-distance.
	idxbase = @spawn Pselection[ AlphaShapes.findClosestPoint(
		P[:, face], P[:, Pselection], metric = "dd"
	) ]
	idxbase=fetch(idxbase)

	# It prevent from adding the same simplex again (cause it has been
	#  determined in a previous recursive call in the stacktrace).
	if idxbase > blackidx
		if DEBUG println("Excluding $face cause simplex already inside.") end
		return nothing
	end

	σ =  sort([face; idxbase])
	if DEBUG @show "Found face" σ end

	# Check the simplex correctness
	radius, center = AlphaShapes.findRadius(P[:, σ], true)

	for i = 1 : size(P, 2)
		if Lar.norm(center - P[:, i]) < radius - 1.e-14
			# @assert i ∉ Pselection "ERROR: Numerical error
			# 	evaluating minimum radius for $σ"
			if DEBUG println("$σ discarded due to a closer point.") end
			return nothing
		end
	end

	return σ
end

#-------------------------------------------------------------------------------

"""
 	firstDeWallSimplex(
		P::Lar.Points, ax::Int64, off::Float64;
		DEBUG = false
	)::Array{Int64,1}

Returns the indices array of the points in `P` that form the first thetrahedron
built over the Wall if the `ax` axes with contant term `off`.
If the keyword argument `DEBUG` is set to true than all the procedure is shown.

# Examples
```jldoctest

julia> V = [
	0.0 1.0 0.0 0.0 2.0
	0.0 0.0 1.0 0.0 2.0
	0.0 0.0 0.0 1.0 2.0
];

julia> firstDeWallSimplex(V, 1, AlphaShapes.findMedian(V,1))
4-element Array{Int64,1}:
 1
 2
 3
 4

```
"""
@timeit to "firstDeWallSimplex" function firstDeWallSimplex(
		P::Lar.Points,
		ax::Int64,
		off::Float64;
		DEBUG = false
	)::Array{Int64,1}

	dim = size(P, 1)
    n = size(P, 2)

	if DEBUG println("Determine first Simplex with ax = $ax") end
    # the first point of the simplex is the one with coordinate `ax` maximal
    #  such that it is less than `off` (closer to α from minus)

	Pselection = findall(x -> x < off, P[ax, :])
	# it gives an error if no point are less than `off`
	#  in fact it means that all the points are located on the median,
	#  with respect to `ax`.
	@assert !isempty(Pselection) "firstDeWallSimplex: not able to build the first Delaunay
		dimplex; all the points have the same `ax` coordinate."

	newidx =  Pselection[findmax(P[ax, Pselection])[2]]
    # indices will store the indices of the simplex ...
    indices = [newidx]                      #Array{Int64,1}
    # ... and `Psimplex` will store the corresponding points
    Psimplex = P[:, newidx][:,:]    #Array{Float64,2}

    # the second point must be seeken across those with coordinate `ax`
    #  grater than `off`
	Pselection = findall(x -> x > off, P[ax, :])

    @inbounds @simd for d = 1 : dim
		idxbase = @spawn AlphaShapes.findClosestPoint(Psimplex, P[:, Pselection])
		idxbase = fetch(idxbase)

		@assert !isnothing(idxbase) "firstDeWallSimplex:
			not able to determine first Delaunay Simplex"
        newidx = Pselection[idxbase]
        indices = [indices; newidx]
        Psimplex = [Psimplex P[:, newidx]]
        Pselection = [i for i = 1 : n if i ∉ indices]
    end

    # Correctness check
	radius, center = AlphaShapes.findRadius(Psimplex, true)

    for i = 1 : n
		@assert Lar.norm(center - P[:, i]) >= radius "firstDeWallSimplex:
			Unable to find first Simplex."
	end

	if DEBUG println("First Simplex = $indices") end

	return indices
end

#-------------------------------------------------------------------------------

"""
	recursiveDelaunayWall(
		P::Lar.Points,
		Pblack::Lar.Points,
		tetraDict::DataStructures.Dict{Array{Int64,1},Array{Float64,1}},
		AFL::Array{Array{Int64,1},1},
		ax::Int64,
		off::Float64,
		positive::Bool;
		DEBUG = false
	)::Lar.Cells

Utility function that prepeares the Divide phase for Delaunay Wall.
Returns the Delaunay Triangulation for the positve or negative subspace of `P`
(according to `positive`) determined by the hyperplane with normal `ax` and
constant term `off`.
If the keyword argument `DEBUG` is set to true than all the procedure is shown.
"""
@timeit to "recursiveDelaunayWall" function recursiveDelaunayWall(
		P::Lar.Points,
		Pblack::Array{Float64},
		tetraDict::DataStructures.Dict{Array{Int64,1},Array{Float64,1}},
		AFL::Array{Array{Int64,1},1},
		ax::Int64,
		off::Float64,
		positive::Bool;
		DEBUG = false
	)::Lar.Cells

	#DEBUG = true

	dim, n = size(P)
	newaxis = mod(ax, dim) + 1

	if DEBUG println("Divide Plus/Minus $positive") end

	Psubset = @spawn findall(x -> (x > off) == positive, P[ax, :])
	Psubset = fetch(Psubset)

	blacklist = @spawn setdiff(unique([(keys(tetraDict)...)...]), Psubset)
	blacklist = fetch(blacklist)

	if !isempty(Pblack)
		Pblack = [Pblack P[:, blacklist]]
	else
	 	Pblack = P[:, blacklist]
	end

	if DEBUG println("Step In") end

	newAFL= @spawn findAFL(Psubset,AFL)
	newAFL=fetch(newAFL)
	newTetraDict= @spawn findTetradict(Psubset,tetraDict)
	newTetraDict=fetch(newTetraDict)
	DT = AlphaShapes.delaunayWall(
			P[:, Psubset],
			newaxis,
			Pblack,
			newAFL,
			newTetraDict,
			DEBUG = DEBUG
		)

	if DEBUG @show "Step Out with " DT end

	return [[Psubset[i] for i in σ] for σ in DT]
end

@timeit to "findTetradict" function findTetradict(Psubset::Array{Int64,1},tetraDict::DataStructures.Dict{Array{Int64,1},Array{Float64,1}})::Dict{Array{Int64,1},Array{Float64,1}}
	return Dict([
		[findall(Psubset.==p)[1] for p in k] => v
		 for (k,v) in tetraDict
			if k ⊆ Psubset
	])
end

@timeit to "findAFL" function findAFL(Psubset::Array{Int64,1},AFL::Array{Array{Int64,1},1})::Array{Array{Int64,1},1}
	return [[findall(Psubset.==p)[1] for p in σ] for σ in AFL]
end

#-------------------------------------------------------------------------------

"""
	updateAFL!(
		P::Lar.Points
		new::Array{Int64,1}[],
		AFLα = Array{Int64,1}[],
		AFLplus = Array{Int64,1}[],
		AFLminus = Array{Int64,1}[],
		ax::Int64, off::Float64;
		DEBUG = false
	)::Bool

Modify the `AFL*` lists of faces by adding to them the faces inside `new`
(that refers to the points `P`) according to their position with respect
to the axis defined by the normal direction `ax` and the contant term `off`.
The function returns a Bool value that states if the operation was succesfully.
If the keyword argument `DEBUG` is set to true than all the procedure is shown.
"""
@timeit to "updateAFL!" function updateAFL!(
		P::Lar.Points,
		newσ::Array{Array{Int64,1},1},
		AFLα::Array{Array{Int64,1},1},
		AFLplus::Array{Array{Int64,1},1},
		AFLminus::Array{Array{Int64,1},1},
		ax::Int64, off::Float64;
		DEBUG = false
	)::Bool

	@simd for face in newσ
		inters = AlphaShapes.planarIntersection(P, face, ax, off)
    	if inters == 0 # intersected by plane α
			AlphaShapes.updatelist!(AFLα, face)
		elseif inters == -1 # in NegHalfspace(α)
        	AlphaShapes.updatelist!(AFLminus, face)
    	elseif inters == 1 # in PosHalfspace(α)
        	AlphaShapes.updatelist!(AFLplus, face)
    	else
			return false
		end
	end

	if DEBUG @show AFLα AFLminus AFLplus end

	return true

end

#-------------------------------------------------------------------------------

"""
	updatelist!(list, element)::Bool

If `element` is in `list` then it is removed (returns `false`);
If `element` is not in `list` then it is added (returns `true`).

# Examples
```jldoctest
julia> list = [[1, 2, 3, 4], [2, 3, 4, 5]]
2-element Array{Array{Int64,1},1}:
 [1, 2, 3, 4]
 [2, 3, 4, 5]

julia> updatelist!(list, [1, 2, 4, 5])
 true

julia> list
3-element Array{Array{Int64,1},1}:
 [1, 2, 3, 4]
 [2, 3, 4, 5]
 [1, 2, 4, 5]

julia> updatelist!(list, [1, 2, 4, 5])
 false

julia> list
2-element Array{Array{Int64,1},1}:
 [1, 2, 3, 4]
 [2, 3, 4, 5]

```
"""

@timeit to "updatelist!" function updatelist!(list, element)::Bool
	if element ∈ list
		setdiff!(list, [element])
		return false
	else
		push!(list, element)
		return true
	end
end

#-------------------------------------------------------------------------------

"""
	updateTetraDict!(
		P::Lar.Points,
		tetraDict::DataStructures.Dict{Array{Int64,1},Array{Float64,1}},
		AFL::Array{Array{Int64,1},1},
		σ::Array{Int64,1}
	)::Nothing

Update the content of `tetraDict` by adding the outer points of the faces of `σ`
in the dictionary.
"""
@timeit to "updateTetraDict" function updateTetraDict!(
		P::Lar.Points,
		tetraDict::DataStructures.Dict{Array{Int64,1},Array{Float64,1}},
		AFL::Array{Array{Int64,1},1},
		σ::Array{Int64,1}
	)::Nothing

	@simd for cell in AFL
		point = setdiff(σ, cell)
		@assert length(point) == 1 "updateTetraDict!: Error during update of TetraDict $σ, $cell"
		tetraDict[ cell ] = P[:, point[1]]
	end

end


text:
fname = joinpath(dirname(@__FILE__), "..", "data", "test.sfp")

s = read_sfp(fname)

@fact length(s)  --> 3
@fact label(s) --> ["Fp1", "Fpz", "Fp2"]
@fact EEG.x(s) --> [-27.747648, -0.085967, 27.676888]
@fact EEG.y(s) --> [98.803864, 103.555275, 99.133354]
@fact EEG.z(s) --> [34.338360, 34.357265, 34.457005]
@fact typeof(s[1].coordinate) --> EEG.Talairach


text:
export mpi_bcast, mpi_init, mpi_recv, mpi_send,
    mpi_sendrecv, mpi_sum,  mpi_finalize, mpi_initialized,
    mpi_finalized, mpi_rank, mpi_size, mpi_sync!, mpi_gather

"""
    mpi_init()

Initialized the MPI session. `mpi_init` must be called before any `run(sess, ...)`.
"""
function mpi_init()
    if mpi_initialized()
        @warn "MPI has already been initialized"
        return
    end
    @eval ccall((:mpi_init, $LIBADCME), Cvoid, ())
end

"""
    mpi_finalize()

Finalize the MPI call.
"""
function mpi_finalize()
    mpi_check()
    @eval ccall((:mpi_finalize, $LIBADCME), Cvoid, ())
end

"""
    mpi_rank()

Returns the rank of current MPI process (rank 0 based).
"""
function mpi_rank()
    mpi_check()
    out = @eval ccall((:mpi_rank, $LIBADCME), Cint, ())
    Int64(out)
end

"""
    mpi_size()

Returns the size of MPI world.
"""
function mpi_size()
    mpi_check()
    @eval ccall((:mpi_size, $LIBADCME), Cint, ())
end

"""
    mpi_finalized()

Returns a boolean indicating whether the current MPI session is finalized.
"""
function mpi_finalized()
    if !mpi_initialized()
        error("MPI has not been initialized.")
    end
    Bool(@eval ccall((:mpi_finalized, $LIBADCME), Cuchar, ()))
end

"""
    mpi_initialized()

Returns a boolean indicating whether the current MPI session is initialized.
"""
function mpi_initialized()
    Bool(@eval ccall((:mpi_initialized, $LIBADCME), Cuchar, ()))
end

"""
    mpi_sync!(message::Array{Int64,1}, root::Int64 = 0)

Sync `message` across all MPI processors.
"""
function mpi_sync!(message::Array{Int64,1}, root::Int64 = 0)
    mpi_check()
    @eval ccall((:mpi_sync, $LIBADCME), Cvoid, (Ptr{Clonglong}, Cint, Cint), $message, Int32(length($message)), Int32($root))
end

function mpi_check()
    if !mpi_initialized()
        error("MPI has not been initialized. Run `mpi_init()` to initialize MPI first.")
    end
    if mpi_finalized()
        error("MPI has been finalized.")
    end
end

function _mpi_sum(a,root::Int64=0)
    mpi_check()
    mpisum_ = load_system_op("mpisum")
    a,root = convert_to_tensor(Any[a,root], [Float64,Int64])
    out = mpisum_(a,root)
    if !isnothing(length(a))
        return set_shape(out, (length(a),))
    else
        return out
    end
end

"""
    mpi_sum(a::Union{Array{Float64}, Float64, PyObject}, root::Int64 = 0)

Sum `a` on the MPI processor `root`.
"""
function mpi_sum(a::Union{Array{Float64}, Float64, PyObject}, root::Int64 = 0)
    a = convert_to_tensor(a, dtype = Float64)
    if length(size(a))==0
        a = reshape(a, (1,))
        out = _mpi_sum(a, root)
        return squeeze(out)
    elseif length(size(a))==1
        return _mpi_sum(a, root)
    elseif nothing in size(a)
        error("The shape of 1st input $(size(a)) contains nothing. mpi_sum is not able to determine
        the output shape.")
    else
        s = size(a)
        a = reshape(a, (-1,))
        out = _mpi_sum(a, root)
        return reshape(out, s)
    end
end

function _mpi_bcast(a,root::Int64=0)
    mpi_check()
    mpibcast_ = load_system_op("mpibcast")
    a,root = convert_to_tensor(Any[a,root], [Float64,Int64])
    out = mpibcast_(a,root)
    if !isnothing(length(a))
        return set_shape(out, (length(a),))
    else
        return out
    end
end

"""
    mpi_bcast(a::Union{Array{Float64}, Float64, PyObject}, root::Int64 = 0)

Broadcast `a` from processor `root` to all other processors.
"""
function mpi_bcast(a::Union{Array{Float64}, Float64, PyObject}, root::Int64 = 0)
    a = convert_to_tensor(a, dtype = Float64)
    if length(size(a))==0
        a = reshape(a, (1,))
        out = _mpi_bcast(a, root)
        return squeeze(out)
    elseif length(size(a))==1
        return _mpi_bcast(a, root)
    elseif nothing in size(a)
        error("The shape of 1st input $(size(a)) contains nothing. mpi_bcast is not able to determine
        the output shape.")
    else
        s = size(a)
        a = reshape(a, (-1,))
        out = _mpi_bcast(a, root)
        return reshape(out, s)
    end
end

function _mpi_send(a,dest::Int64,tag::Int64=0)
    mpisend_ = load_system_op("mpisend")
    a,dest,tag = convert_to_tensor(Any[a,dest,tag], [Float64,Int64,Int64])
    out = mpisend_(a,dest,tag)
    if !isnothing(length(a))
        return set_shape(out, (length(a),))
    else
        return out
    end
end

"""
    mpi_send(a::Union{Array{Float64}, Float64, PyObject}, dest::Int64,root::Int64 = 0)

Sends `a` to processor `dest`. `a` itself is returned so that the send action can be added to the computational graph.
"""
function mpi_send(a::Union{Array{Float64}, Float64, PyObject}, dest::Int64,root::Int64 = 0)
    a = convert_to_tensor(a, dtype = Float64)
    if length(size(a))==0
        a = reshape(a, (1,))
        out = _mpi_send(a,dest, root)
        return squeeze(out)
    elseif length(size(a))==1
        return _mpi_send(a, dest, root)
    elseif nothing in size(a)
        error("The shape of 1st input $(size(a)) contains nothing. mpi_send is not able to determine
        the output shape.")
    else
        s = size(a)
        a = reshape(a, (-1,))
        out = _mpi_send(a, dest, root)
        return reshape(out, s)
    end
end

function _mpi_recv(a,src::Int64,tag::Int64=0)
    mpirecv_ = load_system_op("mpirecv")
    a,src,tag = convert_to_tensor(Any[a,src,tag], [Float64,Int64,Int64])
    mpirecv_(a,src,tag)
end

"""
    mpi_recv(a::Union{Array{Float64}, Float64, PyObject}, src::Int64, tag::Int64 = 0)

Receives an array from processor `src`. `mpi_recv` requires an input for gradient backpropagation.
Typically we can write

```julia
r = mpi_rank()
a = constant(Float64(r))
if r==1
    a = mpi_send(a, 0)
end
if r==0
    a = mpi_recv(a, 1)
end
```

Then `a=1` on both processor 0 and processor 1.
"""
function mpi_recv(a::Union{Array{Float64}, Float64, PyObject}, src::Int64, tag::Int64 = 0)
    a = convert_to_tensor(a, dtype = Float64)
    if length(size(a))==0
        a = reshape(a, (1,))
        out = _mpi_recv(a, src, tag)
        return squeeze(out)
    elseif length(size(a))==1
        return _mpi_recv(a, src, tag)
    elseif nothing in size(a)
        error("The shape of 1st input $(size(a)) contains nothing. mpi_recv is not able to determine
        the output shape.")
    else
        s = size(a)
        a = reshape(a, (-1,))
        out = _mpi_recv(a, src, tag)
        return reshape(out, s)
    end
end

"""
    mpi_sendrecv(a::Union{Array{Float64}, Float64, PyObject}, dest::Int64, src::Int64, tag::Int64=0)

A convenient wrapper for `mpi_send` followed by `mpi_recv`.
"""
function mpi_sendrecv(a::Union{Array{Float64}, Float64, PyObject}, dest::Int64, src::Int64, tag::Int64=0)
    r = mpi_rank()
    @assert src != dest
    if r==src
        a = mpi_send(a, dest, tag)
    elseif r==dest
        a = mpi_recv(a, src, tag)
    end
    a
end

"""
    mpi_gather(u::Union{Array{Float64, 1}, PyObject})

Gathers all the vectors from different processes to the root process. The function returns
a long vector which concatenates of local vectors in the order of process IDs.
"""
function mpi_gather(u::Union{Array{Float64, 1}, PyObject})
    mpigather_ = load_system_op("mpigather")
    u = convert_to_tensor(Any[u], [Float64]); u = u[1]
    out = mpigather_(u)
    set_shape(out, (mpi_size()*length(u)))
end


text:
using PyPlot
using JLD2
function simulate(σ, κ, τ)
    z = zeros(10000)
    α = 0.5
    Δt = 0.01
    z[1] = τ
    for i = 2:length(z)
        x = z[i-1]
        W = randn()
        z[i] = (
            x + κ * (τ - α*x)*Δt + σ * sqrt(x) * sqrt(Δt) * W  + 1/4*σ^2*Δt*(W^2-1)
        )/(1+(1-α)*κ*Δt)
    end
    z
end

close("all")
z = simulate(0.08, 0.5, 0.06)
@save "data.jld2" z
plot(z)
savefig("test.png")


text:
using StructuralCausalModels, Test

ggm_1 = " DAG(U ~X+V, S1~U, W~V+Y, S2~W, order=FALSE)";
ggm_2 = " DAG(U ~X+V, S1~U, W~V+Y, S2~W)";

@testset "ggm conversions" begin

  d1 = from_ggm(ggm_1)
  @test d1 == OrderedDict(
    :U  => [:X, :V],
    :S1 => :U,
    :W  => [:V, :Y],
    :S2 => :W
  )
  d2 = from_ggm(ggm_2)
  @test d2 == OrderedDict(
    :U  => [:X, :V],
    :S1 => :U,
    :W  => [:V, :Y],
    :S2 => :W
  )

  g1 = to_ggm(d1)
  @test g1 == "DAG(U ~ X + V, S1 ~ U, W ~ V + Y, S2 ~ W)"
  g2 = to_ggm(d2)
  @test g2 == "DAG(U ~ X + V, S1 ~ U, W ~ V + Y, S2 ~ W)"

  d3 =OrderedDict(
    [:u, :w] => :v,
    [:u] => :x,
    :s1 => [:u],
    :w => :y,
    [:s2] => [:w]
  );

  g3 = to_ggm(d3)
  @test g3 == "DAG(u ~ v, w ~ v, u ~ x, s1 ~ u, w ~ y, s2 ~ w)"

  g4 = to_ggm(d3, order=true)
  @test g4 == "DAG(u ~ v, w ~ v, u ~ x, s1 ~ u, w ~ y, s2 ~ w, order=TRUE)"

  dag = DAG("ggm_1", ggm_1)
  @test dag.d == OrderedDict(
    :U  => [:X, :V],
    :S1 => :U,
    :W  => [:V, :Y],
    :S2 => :W
  )

end