text:
using BinaryBuilder, Pkg

name = "PROJ"
version = v"7.2.0"

# Collection of sources required to build PROJ
sources = [
    ArchiveSource("https://download.osgeo.org/proj/proj-$version.tar.gz",
        "2957798e5fe295ff96a2af1889d0428e486363d210889422f76dd744f7885763"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/proj-*/

# sqlite needed to build proj.db, so this should not be the
# cross-compiled one since it needs to be executed on the host
apk add sqlite
# Get rid of target sqlite3, to avoid it's picked up by the build system
rm "${bindir}/sqlite3${exeext}"

if [[ ${target} == *mingw* ]]; then
    SQLITE3_LIBRARY=${libdir}/libsqlite3-0.dll
    CURL_LIBRARY=${libdir}/libcurl-4.dll
    TIFF_LIBRARY_RELEASE=${libdir}/libtiff-5.dll
else
    SQLITE3_LIBRARY=${libdir}/libsqlite3.${dlext}
    CURL_LIBRARY=${libdir}/libcurl.${dlext}
    TIFF_LIBRARY_RELEASE=${libdir}/libtiff.${dlext}
fi

mkdir build
cd build
cmake -DCMAKE_INSTALL_PREFIX=$prefix \
      -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
      -DCMAKE_BUILD_TYPE=Release \
      -DBUILD_SHARED_LIBS=ON \
      -DBUILD_TESTING=OFF \
      -DSQLITE3_INCLUDE_DIR=$prefix/include \
      -DSQLITE3_LIBRARY=$SQLITE3_LIBRARY \
      -DCURL_INCLUDE_DIR=$prefix/include \
      -DCURL_LIBRARY=$CURL_LIBRARY \
      -DTIFF_INCLUDE_DIR=$prefix/include \
      -DTIFF_LIBRARY_RELEASE=$TIFF_LIBRARY_RELEASE \
      ..
make -j${nproc}
make install
"""

platforms = expand_cxxstring_abis(supported_platforms())

# The products that we will ensure are always built
products = [
    LibraryProduct(["libproj", "libproj_$(version.major)_$(version.minor)"], :libproj),

    # Excecutables
    ExecutableProduct("cct", :cct),
    ExecutableProduct("cs2cs", :cs2cs),
    ExecutableProduct("geod", :geod),
    ExecutableProduct("gie", :gie),
    ExecutableProduct("proj", :proj),
    ExecutableProduct("projinfo", :projinfo),
    ExecutableProduct("projsync", :projsync),

    # complete contents of share/proj, must be kept up to date
    FileProduct(joinpath("share", "proj", "CH"), :ch),
    FileProduct(joinpath("share", "proj", "GL27"), :gl27),
    FileProduct(joinpath("share", "proj", "ITRF2000"), :itrf2000),
    FileProduct(joinpath("share", "proj", "ITRF2008"), :itrf2008),
    FileProduct(joinpath("share", "proj", "ITRF2014"), :itrf2014),
    FileProduct(joinpath("share", "proj", "nad.lst"), :nad_lst),
    FileProduct(joinpath("share", "proj", "nad27"), :nad27),
    FileProduct(joinpath("share", "proj", "nad83"), :nad83),
    FileProduct(joinpath("share", "proj", "other.extra"), :other_extra),
    FileProduct(joinpath("share", "proj", "proj.db"), :proj_db),
    FileProduct(joinpath("share", "proj", "proj.ini"), :proj_ini),
    FileProduct(joinpath("share", "proj", "world"), :world),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="SQLite_jll", uuid="76ed43ae-9a5d-5a62-8c75-30186b810ce8")),
    Dependency(PackageSpec(name="Libtiff_jll", uuid="89763e89-9b03-5906-acba-b20f662cd828")),
    # libcurl changed compatibility version for macOS from v7.71 to v7.73 (v11
    # to v12)
    Dependency(PackageSpec(name="LibCURL_jll", uuid="deac9b47-8bc7-5906-a0fe-35ac56dc84c0"), v"7.71.1"),
    # The following libraries are dependencies of LibCURL_jll which is now a
    # stdlib, but the stdlib doesn't explicitly list its dependencies
    Dependency("LibSSH2_jll"),
    Dependency("MbedTLS_jll"),
    Dependency("nghttp2_jll"),
    Dependency("Zlib_jll"),
]

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

max_stars_repo_name: vchuravy/Yggdrasil

max_issues_repo_name: vchuravy/Yggdrasil

max_forks_repo_name: vchuravy/Yggdrasil


text:
module Particles

using
    LinearAlgebra,
    Random,
    SpecialFunctions,
    Statistics

using
    Distributions,
    PDMats,
    ConjugatePriors,
    StatsBase,
    StatsFuns

import StatsBase: fit!, fit

using Distributions: NormalStats, MvNormalStats, GenericMvTDist
using ConjugatePriors: posterior_canon, NormalInverseChisq, NormalInverseWishart
using StatsBase: Weights
using StatsFuns: logmvgamma, logπ

export
    ParticleFilter,
    FearnheadParticles,
    ChenLiuParticles,
    InfiniteParticle,
    Component,
    NormalInverseChisq,
    GibbsCRP,
    GibbsCRPSamples,
    ChineseRestaurantProcess,
    StickyCRP,
    ChangePoint,
    NStatePrior,
    Labeled,
    fit,
    fit!,
    putatives,
    weight,
    nobs,
    particles,
    assignments,
    assignment_similarity,
    posterior_predictive,
    marginal_lhood,
    marginal_log_lhood,
    marginal_posterior,
    marginal_log_posterior,
    log_prior,
    marginal_log_prior,
    simulate,
    ncomponents_dist,
    sample!,
    state_entropy,
    randindex

include("component.jl")
include("statepriors.jl")
include("particle.jl")
include("filters.jl")
include("fearnhead.jl")
include("chenliu.jl")
include("gibbs.jl")
include("labeled.jl")

end # module

ext: jl


text:
#=
MOI.Test linear and conic tests
=#

using Test
import MathOptInterface
const MOI = MathOptInterface
const MOIT = MOI.Test
const MOIU = MOI.Utilities
import Hypatia

unit_exclude = [
    "solve_qcp_edge_cases",
    "solve_qp_edge_cases",
    "solve_integer_edge_cases",
    "solve_objbound_edge_cases",
    "solve_zero_one_with_bounds_1",
    "solve_zero_one_with_bounds_2",
    "solve_zero_one_with_bounds_3",
    "solve_unbounded_model", # dual equalities are inconsistent, so detect dual infeasibility but currently no certificate or status
    "number_threads", # no way to set threads currently
    "solve_result_index", # TODO only get one result from Hypatia
    ]

conic_exclude = String[
    # "lin",
    # "sdp",
    # "norminf",
    # "normone",
    # "soc",
    # "rsoc",
    # "normspec",
    # "normnuc",
    # "pow",
    # "dualpow",
    # "geomean",
    # "rootdet",
    "rootdets",
    # "exp",
    # "dualexp",
    # "logdet",
    "logdets",
    # "relentr",
    ]

function test_moi(T::Type{<:Real}; solver_options...)
    optimizer = MOIU.CachingOptimizer(
        MOIU.UniversalFallback(MOIU.Model{T}()),
        Hypatia.Optimizer{T}(; solver_options...)
        )

    tol = 2sqrt(sqrt(eps(T)))
    config = MOIT.TestConfig{T}(
        atol = tol,
        rtol = tol,
        solve = true,
        query = true,
        modify_lhs = true,
        duals = true,
        infeas_certificates = true,
        )

    @testset "linear tests" begin
        MOIT.contlineartest(optimizer, config)
    end

    if T == Float64
        # test other real types, waiting for https://github.com/jump-dev/MathOptInterface.jl/issues/841
        @testset "unit tests" begin
            MOIT.unittest(optimizer, config, unit_exclude)
        end
        @testset "conic tests" begin
            bridged = MOI.Bridges.Constraint.Square{T}(optimizer)
            MOIT.contconictest(bridged, config, conic_exclude)
        end
    end

    return
end


text:
module Kalman

using Distributions, GaussianDistributions
using Trajectories, DynamicIterators
using Random, LinearAlgebra

import DynamicIterators: evolve, dyniterate, @returnnothing
import Base: iterate, IteratorSize, IteratorEltype, eltype, length

import Random.rand
import GaussianDistributions: ⊕

meancov(G) = mean(G), cov(G)
meancov(G::Tuple) = G

symmetrize!(A::Symmetric) = A
symmetrize!(A) = Symmetric(A)
symmetrize!(A::Number) = A

struct Condition{T,S} <: Message2
    u::T
    ll::S
end
struct Filter{T,S} <: Message2
    u::T
    ll::S
end

macro NT(args...)
    :(NamedTuple{($(args)...,)})
end

include("dyniterate.jl")

export LinearEvolution, Evolution
include("evolution.jl")

export LinearObservation, LinearObservationModel, GenericLinearObservationModel, Observation, Observe
include("observation.jl")

export LinearStateSpaceModel, StateObs
include("statespacemodel.jl")

export kalmanfilter
include("filter.jl")

export rts_smoother
include("smoother.jl")
include("backwardsampler.jl")

include("combinator.jl")

end # module


text:
# test the components in optBL separately
# Claudia January 2015

globalerror = false
println("--------- Case G --------------")
include("../examples/case_g_example.jl");

q1 = Quartet(1,["6","7","4","8"],[0.5,0.4,0.1]);
q2 = Quartet(2,["6","7","10","8"],[0.5,0.4,0.1]);
q3 = Quartet(3,["10","7","4","8"],[0.5,0.4,0.1]);
q4 = Quartet(4,["6","10","4","8"],[0.5,0.4,0.1]);
q5 = Quartet(5,["6","7","4","10"],[0.5,0.4,0.1]);

d = DataCF([q1,q2,q3,q4,q5]);
extractQuartet!(net,d)

oldht = net.ht

x = [0.3,1.0,1.5,2.0]
err = false

#println("x is $(x) all changed-----------")

try
    update!(q1.qnet,x,net)
    q1.qnet.edge[3].length !=x[2] || q1.qnet.edge[6].length !=x[3] || q1.qnet.edge[9].length !=x[4] ? error("qnet edges lengths not correct") : nothing
    q1.qnet.edge[5].gamma !=1-x[1] || q1.qnet.edge[7].gamma !=x[1] ? error("qnet edges gammas not correct") : nothing

    update!(q2.qnet,x,net)
    q2.qnet.edge[4].length !=x[3] || q2.qnet.edge[7].length !=x[4] ? error("qnet edges lengths not correct") : nothing
    q2.qnet.edge[3].gamma !=1-x[1] || q2.qnet.edge[5].gamma !=x[1] ? error("qnet edges gammas not correct") : nothing

    update!(q3.qnet,x,net)
    q3.qnet.edge[4].length !=x[3] || q3.qnet.edge[7].length !=x[4] ? error("qnet edges lengths not correct") : nothing
    q3.qnet.edge[3].gamma !=1-x[1] || q3.qnet.edge[5].gamma !=x[1] ? error("qnet edges gammas not correct") : nothing

    update!(q4.qnet,x,net)
    q4.qnet.edge[3].length !=x[2] || q4.qnet.edge[4].length !=x[3] ? error("qnet edges lengths not correct") : nothing

    update!(q5.qnet,x,net)
    q5.qnet.edge[3].length !=x[2] || q5.qnet.edge[6].length !=x[3] ? error("qnet edges lengths not correct") : nothing
    q5.qnet.edge[5].gamma !=1-x[1] || q5.qnet.edge[7].gamma !=x[1] ? error("qnet edges gammas not correct") : nothing

    reduce(&,[q.qnet.changed for q in d.quartet]) || error("all qnet should be changed")

    update!(net,x)
    net.ht == x || ("net.ht not correctly changed to x with update")
catch
    err = true
end

for q in d.quartet
    update!(q.qnet,oldht,net)
end
update!(net,oldht)

x = [0.1,0.2,0.1,2.0] # changing t9 only
#println("x is $(x) changing t9 only-----------")

try
    update!(q1.qnet,x,net)
    q1.qnet.edge[3].length !=x[2] || q1.qnet.edge[6].length !=x[3] || q1.qnet.edge[9].length !=x[4] ? error("qnet edges lengths not correct") : nothing
    q1.qnet.edge[5].gamma !=1-x[1] || q1.qnet.edge[7].gamma !=x[1] ? error("qnet edges gammas not correct") : nothing

    update!(q2.qnet,x,net)
    q2.qnet.edge[4].length !=x[3] || q2.qnet.edge[7].length !=x[4] ? error("qnet edges lengths not correct") : nothing
    q2.qnet.edge[3].gamma !=1-x[1] || q2.qnet.edge[5].gamma !=x[1] ? error("qnet edges gammas not correct") : nothing

    update!(q3.qnet,x,net)
    q3.qnet.edge[4].length !=x[3] || q3.qnet.edge[7].length !=x[4] ? error("qnet edges lengths not correct") : nothing
    q3.qnet.edge[3].gamma !=1-x[1] || q3.qnet.edge[5].gamma !=x[1] ? error("qnet edges gammas not correct") : nothing

    update!(q4.qnet,x,net)
    q4.qnet.edge[3].length !=x[2] || q4.qnet.edge[4].length !=x[3] ? error("qnet edges lengths not correct") : nothing

    update!(q5.qnet,x,net)
    q5.qnet.edge[3].length !=x[2] || q5.qnet.edge[6].length !=x[3] ? error("qnet edges lengths not correct") : nothing
    q5.qnet.edge[5].gamma !=1-x[1] || q5.qnet.edge[7].gamma !=x[1] ? error("qnet edges gammas not correct") : nothing

    [q.qnet.changed for q in d.quartet] == [true,true,true,false,false] || error("q.qnet.changed not correct for all quartets")

    update!(net,x)
    net.ht == x || error("net.ht not correctly changed to x with update")
catch
    err = true
end

for q in d.quartet
    update!(q.qnet,oldht,net)
end
update!(net,oldht)

x = [0.3,0.2,0.1,2.0] # changing gamma and t9 only
#println("x is $(x) changing gamma and t9 only-----------")

try
    update!(q1.qnet,x,net)
    q1.qnet.edge[3].length !=x[2] || q1.qnet.edge[6].length !=x[3] || q1.qnet.edge[9].length !=x[4] ? error("qnet edges lengths not correct") : nothing
    q1.qnet.edge[5].gamma !=1-x[1] || q1.qnet.edge[7].gamma !=x[1] ? error("qnet edges gammas not correct") : nothing

    update!(q2.qnet,x,net)
    q2.qnet.edge[4].length !=x[3] || q2.qnet.edge[7].length !=x[4] ? error("qnet edges lengths not correct") : nothing
    q2.qnet.edge[3].gamma !=1-x[1] || q2.qnet.edge[5].gamma !=x[1] ? error("qnet edges gammas not correct") : nothing

    update!(q3.qnet,x,net)
    q3.qnet.edge[4].length !=x[3] || q3.qnet.edge[7].length !=x[4] ? error("qnet edges lengths not correct") : nothing
    q3.qnet.edge[3].gamma !=1-x[1] || q3.qnet.edge[5].gamma !=x[1] ? error("qnet edges gammas not correct") : nothing

    update!(q4.qnet,x,net)
    q4.qnet.edge[3].length !=x[2] || q4.qnet.edge[4].length !=x[3] ? error("qnet edges lengths not correct") : nothing

    update!(q5.qnet,x,net)
    q5.qnet.edge[3].length !=x[2] || q5.qnet.edge[6].length !=x[3] ? error("qnet edges lengths not correct") : nothing
    q5.qnet.edge[5].gamma !=1-x[1] || q5.qnet.edge[7].gamma !=x[1] ? error("qnet edges gammas not correct") : nothing

    [q.qnet.changed for q in d.quartet] == [true,true,true,false,true] || error("q.qnet.changed not correct for all quartets")

    update!(net,x)
    net.ht == x || error("net.ht not correctly changed to x with update")
catch
    err = true
end

for q in d.quartet
    update!(q.qnet,oldht,net)
end
update!(net,oldht)

# ---- calculateExpCF
x = [0.3,0.2,0.1,2.0] # changing gamma and t9 only
#println("---- calculate expCF for $(x)")
try
    calculateExpCFAll!(d,x,net)
    reduce(&,map(approxEq,q1.qnet.expCF,[(1-x[1])/3*exp(-x[2])+x[1]/3*exp(-x[2]-x[3]-x[4]),(1-x[1])*(1-2/3*exp(-x[2]))+x[1]*(1-2/3*exp(-x[2]-x[3]-x[4])),
                                    (1-x[1])/3*exp(-x[2])+x[1]/3*exp(-x[2]-x[3]-x[4])])) || error("q1 expCF wrong")
    reduce(&,map(approxEq,q2.qnet.expCF, [(1-x[1])*(1-2/3*exp(-x[3]))+x[1]*(1/3*exp(-x[4])),(1-x[1])*(1/3*exp(-x[3]))+x[1]*(1-2/3*exp(-x[4])),
                                     (1-x[1])/3*exp(-x[3])+x[1]/3*exp(-x[4])])) || error("q2 expCF wrong")
    reduce(&,map(approxEq,q3.qnet.expCF, [(1-x[1])/3*exp(-x[3])+x[1]/3*exp(-x[4]),(1-x[1])*(1/3*exp(-x[3]))+x[1]*(1-2/3*exp(-x[4])),(1-x[1])*(1-2/3*exp(-x[3]))+x[1]*(1/3*exp(-x[4]))])) || error("q3 expCF wrong")
    reduce(&,map(approxEq,q4.qnet.expCF,[1/3*exp(-x[2]-x[3]),1-2/3*exp(-x[2]-x[3]),1/3*exp(-x[2]-x[3])])) || error("q4 expCF wrong")
    reduce(&,map(approxEq,q5.qnet.expCF,[(1-x[1])/3*exp(-x[2])+x[1]/3*exp(-x[2]-x[3]),(1-x[1])*(1-2/3*exp(-x[2]))+x[1]*(1-2/3*exp(-x[2]-x[3])),
                                    (1-x[1])/3*exp(-x[2])+x[1]/3*exp(-x[2]-x[3])])) || error("q5 expCF wrong")
catch
    err = true
end

#logPseudoLik(d)

if !err
    #println("-------------Case G: NO ERRORS!------------")
else
    globalerror = true
end

println("--------- Case F Bad Diamond I --------------")
include("../examples/case_f_example.jl");
parameters!(net)

q1 = Quartet(1,["6","7","4","8"],[0.5,0.4,0.1]);
q2 = Quartet(2,["6","7","10","8"],[0.5,0.4,0.1]);
q3 = Quartet(3,["10","7","4","8"],[0.5,0.4,0.1]);
q4 = Quartet(4,["6","10","4","8"],[0.5,0.4,0.1]);
q5 = Quartet(5,["6","7","4","10"],[0.5,0.4,0.1]);

d = DataCF([q1,q2,q3,q4,q5]);
extractQuartet!(net,d)

oldht = net.ht

x = [0.4,0.2,0.1]
err = false

#println("x is $(x) all changed-----------")

try
    update!(q1.qnet,x,net)
    q1.qnet.node[1].gammaz !=x[2] || q1.qnet.node[3].gammaz !=x[3]  ? error("qnet gammaz not correct") : nothing

    update!(q2.qnet,x,net)
    q2.qnet.edge[4].length !=x[1]  ? error("qnet edges lengths not correct") : nothing

    update!(q3.qnet,x,net)
    q3.qnet.edge[5].length !=x[1] ? error("qnet edges lengths not correct") : nothing
    q3.qnet.edge[3].length !=-log(1-x[3]) ? error("qnet edges gammaz not correct") : nothing

    update!(q4.qnet,x,net)
    q4.qnet.edge[5].length !=x[1] ? error("qnet edges lengths not correct") : nothing
    q4.qnet.edge[3].length !=-log(1-x[2]) ? error("qnet edges gammaz not correct") : nothing

    update!(q5.qnet,x,net)
    q1.qnet.node[1].gammaz !=x[2] || q1.qnet.node[3].gammaz !=x[3]  ? error("qnet gammaz not correct") : nothing

    reduce(&,[q.qnet.changed for q in d.quartet]) || error("all qnet should be changed")

    update!(net,x)
    net.ht == x || ("net.ht not correctly changed to x with update")
catch
    err = true
end

for q in d.quartet
    update!(q.qnet,oldht,net)
end
update!(net,oldht)

x = [0.1,0.2,0.1] # changing gammaz1, gammaz2 only
#println("x is $(x) changing gammaz1, gammaz2 only-----------")

try
    update!(q1.qnet,x,net)
    q1.qnet.node[1].gammaz !=x[2] || q1.qnet.node[3].gammaz !=x[3]  ? error("qnet gammaz not correct") : nothing

    update!(q2.qnet,x,net)
    q2.qnet.edge[4].length !=x[1]  ? error("qnet edges lengths not correct") : nothing

    update!(q3.qnet,x,net)
    q3.qnet.edge[5].length !=x[1] ? error("qnet edges lengths not correct") : nothing
    q3.qnet.edge[3].length !=-log(1-x[3]) ? error("qnet edges gammaz not correct") : nothing

    update!(q4.qnet,x,net)
    q4.qnet.edge[5].length !=x[1] ? error("qnet edges lengths not correct") : nothing
    q4.qnet.edge[3].length !=-log(1-x[2]) ? error("qnet edges gammaz not correct") : nothing

    update!(q5.qnet,x,net)
    q1.qnet.node[1].gammaz !=x[2] || q1.qnet.node[3].gammaz !=x[3]  ? error("qnet gammaz not correct") : nothing

    [q.qnet.changed for q in d.quartet] == [true,false,true,true,true] || error("q.qnet.changed not correct for all quartets")

    update!(net,x)
    net.ht == x || error("net.ht not correctly changed to x with update")
catch
    err = true
end

for q in d.quartet
    update!(q.qnet,oldht,net)
end
update!(net,oldht)

# ---- calculateExpCF
x = [0.3,0.2,0.1]
#println("---- calculate expCF for $(x)")
try
    calculateExpCFAll!(d,x,net)
    reduce(&,map(approxEq,q1.qnet.expCF,[(1-x[2]-x[3])/3,(1+2*x[2]-x[3])/3,(1-x[2]+2*x[3])/3])) || error("q1 expCF wrong")
    reduce(&,map(approxEq,q2.qnet.expCF, [1-2/3*exp(-x[1]),1/3*exp(-x[1]),1/3*exp(-x[1])])) || error("q2 expCF wrong")
    reduce(&,map(approxEq,q3.qnet.expCF, [1/3*exp(-x[1]+log(1-x[3])),1/3*exp(-x[1]+log(1-x[3])),1-2/3*exp(-x[1]+log(1-x[3]))])) || error("q3 expCF wrong")
    reduce(&,map(approxEq,q4.qnet.expCF,[1/3*exp(-x[1]+log(1-x[2])),1-2/3*exp(-x[1]+log(1-x[2])),1/3*exp(-x[1]+log(1-x[2]))])) || error("q4 expCF wrong")
    reduce(&,map(approxEq,q5.qnet.expCF,[(1-x[2]-x[3])/3,(1+2*x[2]-x[3])/3,(1-x[2]+2*x[3])/3])) || error("q5 expCF wrong")

catch
    err = true
end

#logPseudoLik(d)

if !err
    #println("-------------Case F: NO ERRORS!------------")
else
    globalerror = true
end

println("--------- Case I Bad Diamond II --------------")
include("../examples/case_i_example.jl");

q1 = Quartet(1,["6","7","4","8"],[0.5,0.4,0.1]);
q2 = Quartet(2,["6","7","10","8"],[0.5,0.4,0.1]);
q3 = Quartet(3,["10","7","4","8"],[0.5,0.4,0.1]);
q4 = Quartet(4,["6","10","4","8"],[0.5,0.4,0.1]);
q5 = Quartet(5,["6","7","4","10"],[0.5,0.4,0.1]);

d = DataCF([q1,q2,q3,q4,q5]);
extractQuartet!(net,d)

oldht = net.ht

x = [0.2,0.2,0.1,0.1,0.1]
err = false

#println("x is $(x) all changed-----------")

try
    update!(q1.qnet,x,net)
    q1.qnet.edge[7].gamma != x[1] || q1.qnet.edge[2].gamma != 1-x[1]  ? error("qnet gamma not correct") : nothing
    q1.qnet.edge[4].length != x[3] || q1.qnet.edge[8].length != x[4] ? error("qnet edge lengths not correct") : nothing

    update!(q2.qnet,x,net)
    q2.qnet.edge[8].gamma != x[1] || q2.qnet.edge[4].gamma != 1-x[1]  ? error("qnet gamma not correct") : nothing
    (q2.qnet.edge[4].length !=x[2] || q2.qnet.edge[6].length !=x[3] || q2.qnet.edge[8].length !=x[4] || q2.qnet.edge[9].length !=x[5])  ? error("qnet edges lengths not correct") : nothing

    update!(q3.qnet,x,net)
    q3.qnet.edge[8].gamma != x[1] || q3.qnet.edge[4].gamma != 1-x[1]  ? error("qnet gamma not correct") : nothing
    (q2.qnet.edge[4].length !=x[2] || q2.qnet.edge[6].length !=x[3] || q2.qnet.edge[8].length !=x[4] || q2.qnet.edge[9].length !=x[5])  ? error("qnet edges lengths not correct") : nothing

    update!(q4.qnet,x,net)
    q4.qnet.edge[8].gamma != x[1] || q4.qnet.edge[4].gamma != 1-x[1]  ? error("qnet gamma not correct") : nothing
    (q2.qnet.edge[4].length !=x[2] || q2.qnet.edge[6].length !=x[3] || q2.qnet.edge[8].length !=x[4] || q2.qnet.edge[9].length !=x[5])  ? error("qnet edges lengths not correct") : nothing

    update!(q5.qnet,x,net)
    q5.qnet.edge[7].gamma != x[1] || q5.qnet.edge[2].gamma != 1-x[1]  ? error("qnet gamma not correct") : nothing
    q5.qnet.edge[4].length != x[3] || q5.qnet.edge[8].length != x[4] ? error("qnet edge lengths not correct") : nothing

    reduce(&,[q.qnet.changed for q in d.quartet]) || error("all qnet should be changed")

    update!(net,x)
    net.ht == x || ("net.ht not correctly changed to x with update")
catch
    err = true
end

for q in d.quartet
    update!(q.qnet,oldht,net)
end
update!(net,oldht)

x = [0.1,0.2,1.,1.,1.] # changing t3 only
err = false

#println("x is $(x) changing t3 only-----------")

try
    update!(q1.qnet,x,net)
    q1.qnet.edge[7].gamma != x[1] || q1.qnet.edge[2].gamma != 1-x[1]  ? error("qnet gamma not correct") : nothing
    q1.qnet.edge[4].length != x[3] || q1.qnet.edge[8].length != x[4] ? error("qnet edge lengths not correct") : nothing

    update!(q2.qnet,x,net)
    q2.qnet.edge[8].gamma != x[1] || q2.qnet.edge[4].gamma != 1-x[1]  ? error("qnet gamma not correct") : nothing
    (q2.qnet.edge[4].length !=x[2] || q2.qnet.edge[6].length !=x[3] || q2.qnet.edge[8].length !=x[4] || q2.qnet.edge[9].length !=x[5])  ? error("qnet edges lengths not correct") : nothing

    update!(q3.qnet,x,net)
    q3.qnet.edge[8].gamma != x[1] || q3.qnet.edge[4].gamma != 1-x[1]  ? error("qnet gamma not correct") : nothing
    (q2.qnet.edge[4].length !=x[2] || q2.qnet.edge[6].length !=x[3] || q2.qnet.edge[8].length !=x[4] || q2.qnet.edge[9].length !=x[5])  ? error("qnet edges lengths not correct") : nothing

    update!(q4.qnet,x,net)
    q4.qnet.edge[8].gamma != x[1] || q4.qnet.edge[4].gamma != 1-x[1]  ? error("qnet gamma not correct") : nothing
    (q2.qnet.edge[4].length !=x[2] || q2.qnet.edge[6].length !=x[3] || q2.qnet.edge[8].length !=x[4] || q2.qnet.edge[9].length !=x[5])  ? error("qnet edges lengths not correct") : nothing

    update!(q5.qnet,x,net)
    q5.qnet.edge[7].gamma != x[1] || q5.qnet.edge[2].gamma != 1-x[1]  ? error("qnet gamma not correct") : nothing
    q5.qnet.edge[4].length != x[3] || q5.qnet.edge[8].length != x[4] ? error("qnet edge lengths not correct") : nothing

    [q.qnet.changed for q in d.quartet] == [false, true, true, true, false] || error("not all qnet should be changed")

    update!(net,x)
    net.ht == x || ("net.ht not correctly changed to x with update")
catch
    err = true
end

for q in d.quartet
    update!(q.qnet,oldht,net)
end
update!(net,oldht)

# ---- calculateExpCF
x = [0.2,0.2,0.1,0.1,0.1]
#println("---- calculate expCF for $(x)")
try
    calculateExpCFAll!(d,x,net)
    reduce(&,map(approxEq,q1.qnet.expCF,[(1-x[1])*(1/3*exp(-x[3]))+x[1]*(1-2/3*exp(-x[4])),(1-x[1])*(1-2/3*exp(-x[3]))+x[1]*(1/3*exp(-x[4])),(1-x[1])*(1/3*exp(-x[3]))+x[1]*(1/3*exp(-x[4]))])) || error("q1 expCF wrong")
    t=-log(1+x[1]*(1-exp(-x[3]))-x[1]*x[1]*(1-exp(-x[5]-x[4]))-x[1]*x[1]*(1-exp(-x[3]))-(1-x[1])*(1-x[1])*(1-exp(-x[2])))
    reduce(&,map(approxEq,q2.qnet.expCF, [1-2/3*exp(-t),1/3*exp(-t),1/3*exp(-t)])) || error("q2 expCF wrong")
    t=-log(1+x[1]*(1-exp(-x[3]-x[5]))-x[1]*x[1]*(1-exp(-x[4]))-x[1]*x[1]*(1-exp(-x[3]-x[5]))-(1-x[1])*(1-x[1])*(1-exp(-x[2])))
    reduce(&,map(approxEq,q3.qnet.expCF, [1/3*exp(-t),1/3*exp(-t),1-2/3*exp(-t)])) || error("q3 expCF wrong")
    t=-log(1+x[1]*(1-exp(-x[5]))-x[1]*x[1]*(1-exp(-x[5]))-x[1]*x[1]*(1-exp(-x[4]))-(1-x[1])*(1-x[1])*(1-exp(-x[3]-x[2])))
    reduce(&,map(approxEq,q4.qnet.expCF, [1/3*exp(-t),1-2/3*exp(-t),1/3*exp(-t)])) || error("q4 expCF wrong")
    reduce(&,map(approxEq,q5.qnet.expCF,[(1-x[1])*(1/3*exp(-x[3]))+x[1]*(1-2/3*exp(-x[4])),(1-x[1])*(1-2/3*exp(-x[3]))+x[1]*(1/3*exp(-x[4])),(1-x[1])*(1/3*exp(-x[3]))+x[1]*(1/3*exp(-x[4]))])) || error("q5 expCF wrong")

catch
    err = true
end

#logPseudoLik(d)

if !err
    #println("-------------Case I: NO ERRORS!------------")
else
    globalerror = true
end

@test !globalerror


text:
mutable struct ZWarpedArray{TO,TI,N,NC,AA<:AbstractArray} <: AbstractCachedArray{TO,TI,N,NC,AA}
    parent::AA
    tfms::Vector{CoordinateTransformations.Transformation}
    cached::Array{TO,NC}
    current_I::Tuple
    correct_bias::Bool
    sqrt_tfm::Bool
end

cache(A::ZWarpedArray) = A.cached
current_I(A::ZWarpedArray) = A.current_I
set_I!(A::ZWarpedArray, inds) = A.current_I = inds

function ZWarpedArray(img::AbstractArray{T, N}, tfms, out_type=Float64, nd_cache::Int=2; correct_bias=true, sqrt_tfm=false) where {T, N}
    if !(2 < ndims(img) < 5)
        error("Only 3D and 4D images are supported")
    end
    if size(img,3) !== length(tfms)
        error("Input image size in the Z-slice dimension (3) does not match the number of transforms provided")
    end
    if nd_cache < 2 || nd_cache > 3
        error("nd_cache must be 2 or 3")
    end
    if sqrt_tfm && !correct_bias
        warn("Square root transformation will yield incorrect results if bias was not subtracted first")
    end
    if correct_bias && !in(T, (Normed{UInt16, 16}, UInt16))
        warn("Unable to guess the camera bias value for element type $T.  Defaulting to 100/(2^16)")
    end
    cache_rngs = axes(img)[nd_cache+1:ndims(img)]
    ci = map(first, cache_rngs)
    cached = similar(view(img, axes(img)[1:nd_cache]..., ci...), out_type)
    za = ZWarpedArray{out_type, T, ndims(img), nd_cache, typeof(img)}(img, tfms, cached, (ci...,), correct_bias, sqrt_tfm)
    update_cache!(za, (ci...,))
    return za
end

function _update_cache!(c::AbstractArray{T,2}, parent, tfm::CoordinateTransformations.Transformation, pp::Function) where {T}
    c .= pp.(warp_and_resample(parent, tfm)) #TODO: do this in-place
    return nothing
end

function _update_cache!(c::AbstractArray{T,3}, parent, tfms, pp::Function) where {T}
    for (i,v) in enumerate(axes(c,3))
        _update_cache!(view(c, :, :, v), view(parent, :, :, v), tfms[i], pp)
    end
    return nothing
end

function biasval(t::Type)
    native = UInt16(100)
    return ifelse(t==UInt16, native, reinterpret(Normed{UInt16, 16}, native))
end

function update_cache!(A::ZWarpedArray{TO,TI,N,NC,AA}, inds::NTuple{N2,Int}) where {TO,TI,N,NC,AA,N2}
    pslice = view(parent(A), cached_axes(A)..., inds...) #TODO: may be faster to use copy instead of view here
    pp = identity
    #note: the max statements below only seem necessary due to the way fixed point numbers are converted,
    #see Interpolations #118
    if A.correct_bias && A.sqrt_tfm
        pp = x-> sqrt(Float64(max(0, x-biasval(TI))))
    elseif A.correct_bias
        pp = x-> Float64(max(0, x- biasval(TI)))
    elseif A.sqrt_tfm
        pp = x-> sqrt(Float64(x))
    end
    if length(cached_axes(A)) == 2
        tfmidx = findfirst(isequal(first(inds)), axes(A,3))
        _update_cache!(cache(A), pslice, A.tfms[tfmidx], pp)
    else
        _update_cache!(cache(A), pslice, A.tfms, pp)
    end
    set_I!(A, inds)
    return nothing
end

show(io::IO, A::ZWarpedArray{TO}) where {TO} = print(io, "ZWarpedArray of size $(size(A)) mapped to element type $TO\n")
show(io::IO, ::MIME"text/plain", A::ZWarpedArray{TO}) where {TO} = show(io, A)

ZWarpedArray(img::ImageMeta, tfms, out_type=Float64; kwargs...) = ImageMeta(ZWarpedArray(data(img), tfms, out_type; kwargs...), properties(img))
ZWarpedArray(img::AxisArray, tfms, out_type=Float64; kwargs...) = match_axisspacing(ZWarpedArray(data(img),tfms,out_type; kwargs...), img)


text:
using motion_planning
using Test

println("Testing...")


text:
@testset "shah" begin
    img = generate_surface(SynthSphere(50), 0.5)
    Z = Shah()(img)
    @test Z[76,76] ≈ 1.1326467103582539e20
    @test Z[1,1] ≈ 1.1483728689832028e20
    @test Z[151,151] ≈ 1.1450402052589486e20
    @test isapprox(Z[90,90],  1.038114587065121e20, atol = 1e18)
    @test Z[101,36] ≈ 6.890057069914122e19
    @test Z[22, 136] ≈ 1.147427113061455e20
    @test Z[42,42] ≈ 7.171417157958961e19
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiAuthorizationV1beta1SubjectAccessReviewStatus <: SwaggerModel
    allowed::Any # spec type: Union{ Nothing, Bool } # spec name: allowed
    denied::Any # spec type: Union{ Nothing, Bool } # spec name: denied
    evaluationError::Any # spec type: Union{ Nothing, String } # spec name: evaluationError
    reason::Any # spec type: Union{ Nothing, String } # spec name: reason

    function IoK8sApiAuthorizationV1beta1SubjectAccessReviewStatus(;allowed=nothing, denied=nothing, evaluationError=nothing, reason=nothing)
        o = new()
        validate_property(IoK8sApiAuthorizationV1beta1SubjectAccessReviewStatus, Symbol("allowed"), allowed)
        setfield!(o, Symbol("allowed"), allowed)
        validate_property(IoK8sApiAuthorizationV1beta1SubjectAccessReviewStatus, Symbol("denied"), denied)
        setfield!(o, Symbol("denied"), denied)
        validate_property(IoK8sApiAuthorizationV1beta1SubjectAccessReviewStatus, Symbol("evaluationError"), evaluationError)
        setfield!(o, Symbol("evaluationError"), evaluationError)
        validate_property(IoK8sApiAuthorizationV1beta1SubjectAccessReviewStatus, Symbol("reason"), reason)
        setfield!(o, Symbol("reason"), reason)
        o
    end
end # type IoK8sApiAuthorizationV1beta1SubjectAccessReviewStatus

const _property_map_IoK8sApiAuthorizationV1beta1SubjectAccessReviewStatus = Dict{Symbol,Symbol}(Symbol("allowed")=>Symbol("allowed"), Symbol("denied")=>Symbol("denied"), Symbol("evaluationError")=>Symbol("evaluationError"), Symbol("reason")=>Symbol("reason"))
const _property_types_IoK8sApiAuthorizationV1beta1SubjectAccessReviewStatus = Dict{Symbol,String}(Symbol("allowed")=>"Bool", Symbol("denied")=>"Bool", Symbol("evaluationError")=>"String", Symbol("reason")=>"String")
Base.propertynames(::Type{ IoK8sApiAuthorizationV1beta1SubjectAccessReviewStatus }) = collect(keys(_property_map_IoK8sApiAuthorizationV1beta1SubjectAccessReviewStatus))
Swagger.property_type(::Type{ IoK8sApiAuthorizationV1beta1SubjectAccessReviewStatus }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApiAuthorizationV1beta1SubjectAccessReviewStatus[name]))}
Swagger.field_name(::Type{ IoK8sApiAuthorizationV1beta1SubjectAccessReviewStatus }, property_name::Symbol) =  _property_map_IoK8sApiAuthorizationV1beta1SubjectAccessReviewStatus[property_name]

function check_required(o::IoK8sApiAuthorizationV1beta1SubjectAccessReviewStatus)
    (getproperty(o, Symbol("allowed")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiAuthorizationV1beta1SubjectAccessReviewStatus }, name::Symbol, val)
end


text:
@inline function mulsizeof(::Type{T}, x::Number) where {T}
  (Base.allocatedinline(T) ? sizeof(T) : sizeof(Int)) * x
end
@generated function mulsizeof(::Type{T}, ::StaticInt{N}) where {T,N}
  st = Base.allocatedinline(T) ? sizeof(T) : sizeof(Int)
  Expr(:block, Expr(:meta,:inline), Expr(:call, Expr(:curly, :StaticInt, N*st)))
end
@inline mulsizeof(::Type{T}, ::Tuple{}) where {T} = ()
@inline mulsizeof(::Type{T}, x::Tuple{X}) where {T,X} = (mulsizeof(T, getfield(x,1,false)), )
@inline mulsizeof(::Type{T}, x::Tuple{X1,X2,Vararg}) where {T,X1,X2} = (mulsizeof(T, getfield(x,1,false)), mulsizeof(T, Base.tail(x))...)

@inline bytestrides(A::AbstractArray{T}) where {T} = mulsizeof(T, ArrayInterface.strides(A))

@inline memory_reference(A::AbstractArray) = memory_reference(ArrayInterface.device(A), A)
@inline memory_reference(A::BitArray) = Base.unsafe_convert(Ptr{Bit}, A.chunks), A.chunks
@inline memory_reference(::CPUPointer, A) = pointer(A), preserve_buffer(A)
@inline memory_reference(::CPUPointer, A::Union{LinearAlgebra.Adjoint, Base.ReshapedArray, Base.PermutedDimsArray, LinearAlgebra.Transpose}) = memory_reference(CPUPointer(), parent(A))
@inline function memory_reference(::CPUPointer, A::Base.ReinterpretArray{T}) where {T}
  p, m = memory_reference(CPUPointer(), parent(A))
  reinterpret(Ptr{T}, p), m
end
@inline ind_diff(::Base.Slice, ::Any) = Zero()
@inline ind_diff(x::AbstractRange, o) = static_first(x) - o
@inline ind_diff(x::Integer, o) = x - o
@inline memory_reference(::CPUPointer, A::SubArray) = memory_reference_subarray(CPUPointer(), A)
@inline memory_reference(::CPUTuple, A::SubArray) = memory_reference_subarray(CPUTuple(), A)
@inline function memory_reference_subarray(::PT, A::SubArray) where {PT}
  p, m = memory_reference(PT(), parent(A))
  pA = parent(A)
  offset = ArrayInterface.reduce_tup(+, _map(*, _map(ind_diff, A.indices, offsets(pA)), strides(pA)))
  p + sizeof(eltype(A))*offset, m
end
@inline function memory_reference(::CPUTuple, A)
  r = Ref(A)
  Base.unsafe_convert(Ptr{eltype(A)}, Base.pointer_from_objref(r)), r
end
@inline function memory_reference(::ArrayInterface.CheckParent, A)
  P = parent(A)
  if P === A
    memory_reference(ArrayInterface.CPUIndex(), A)
  else
    memory_reference(ArrayInterface.device(P), P)
  end
end
@inline memory_reference(::ArrayInterface.CPUIndex, A) = throw("Memory access for $(typeof(A)) not implemented yet.")

@inline ArrayInterface.contiguous_axis(::Type{A}) where {T,N,C,A<:AbstractStridedPointer{T,N,C}} = StaticInt{C}()
@inline ArrayInterface.contiguous_batch_size(::Type{A}) where {T,N,C,B,A<:AbstractStridedPointer{T,N,C,B}} = StaticInt{B}()
@inline ArrayInterface.stride_rank(::Type{A}) where {T,N,C,B,R,A<:AbstractStridedPointer{T,N,C,B,R}} = _map(StaticInt, R)
@inline memory_reference(A::AbstractStridedPointer) = pointer(A), nothing

@inline Base.eltype(::AbstractStridedPointer{T}) where {T} = T

struct StridedPointer{T,N,C,B,R,X,O} <: AbstractStridedPointer{T,N,C,B,R,X,O}
  p::Ptr{T}
  si::StrideIndex{N,R,C,X,O}
end
@inline stridedpointer(p::Ptr{T}, si::StrideIndex{N,R,C,X,O}, ::StaticInt{B}) where {T,N,R,C,B,X,O} = StridedPointer{T,N,C,B,R,X,O}(p, si)
@inline stridedpointer(p::Ptr{T}, si::StrideIndex{N,R,C,X,O}) where {T,N,R,C,X,O} = StridedPointer{T,N,C,0,R,X,O}(p, si)

@inline bytestrideindex(A::AbstractArray{T}) where {T} = mulstrides(T, StrideIndex(A))

@inline function mulstrides(::Type{T}, si::StrideIndex{N,R,C,X,O}) where {N,R,C,X,O,T}
  StrideIndex{N,R,C}(mulsizeof(T, si.strides), si.offsets)
end
@inline function stridedpointer(A::AbstractArray)
  p, r = memory_reference(A)
  stridedpointer(p, bytestrideindex(A), ArrayInterface.contiguous_batch_size(A))
end
@inline function stridedpointer_preserve(A::AbstractArray)
  p, r = memory_reference(A)

  stridedpointer(p, bytestrideindex(A), ArrayInterface.contiguous_batch_size(A)), r
end
@inline val_stride_rank(::AbstractStridedPointer{T,N,C,B,R}) where {T,N,C,B,R} = Val{R}()
@generated val_dense_dims(::AbstractStridedPointer{T,N}) where {T,N} = Val{ntuple(==(0), Val(N))}()
@inline val_stride_rank(A) = Val(known(stride_rank(A)))
@inline val_dense_dims(A) = Val(known(ArrayInterface.dense_dims(A)))

function zerotupleexpr(N::Int)
  t = Expr(:tuple);
  for n in 1:N
    push!(t.args, :(Zero()))
  end
  Expr(:block, Expr(:meta,:inline), t)
end
@generated zerotuple(::Val{N}) where {N} = zerotupleexpr(N)
# @inline center(p::AbstractStridedPointer{T,N}) where {T,N} = gesp(p, zerotuple(Val(N)))
@inline zero_offsets(si::StrideIndex{N,R,C}) where {N,R,C} = StrideIndex{N,R,C}(getfield(si,:strides), zerotuple(Val(N)))
@inline zero_offsets(sptr::AbstractStridedPointer) = stridedpointer(pointer(sptr), zero_offsets(StrideIndex(sptr)), contiguous_batch_size(sptr))
@inline zstridedpointer(A) = zero_offsets(stridedpointer(A))
@inline function zstridedpointer_preserve(A::AbstractArray{T,N}) where {T,N}
  strd = mulsizeof(T, ArrayInterface.strides(A))
  si = StrideIndex{N,known(stride_rank(A)),Int(contiguous_axis(A))}(strd, zerotuple(Val(N)), Zero())
  p, r = memory_reference(A)
  stridedpointer(p, si, contiguous_batch_size(A)), r
end

@inline Base.similar(sptr::StridedPointer, ptr::Ptr) = stridedpointer(ptr, StrideIndex(sptr), contiguous_batch_size(sptr))

Base.unsafe_convert(::Type{Ptr{T}}, ptr::AbstractStridedPointer{T}) where {T} = pointer(ptr)
# Shouldn't need to special case Array
@generated function nopromote_axis_indicator(::AbstractStridedPointer{<:Any,N}) where {N}
  t = Expr(:tuple); foreach(n -> push!(t.args, True()), 1:N)
  Expr(:block, Expr(:meta, :inline), t)
end

@inline dynamic_offsets(si::StrideIndex{N,R,C}) where {N,R,C} = StrideIndex{N,R,C}(strides(si), _map(Int, offsets(si)))
struct StridedBitPointer{N,C,B,R,X,O} <: AbstractStridedPointer{Bit,N,C,B,R,X,O}
  p::Ptr{Bit}
  si::StrideIndex{N,R,C,X,O}
end
@inline stridedpointer(p::Ptr{Bit}, si::StrideIndex{N,R,C,X,O}, ::StaticInt{B}) where {N,R,C,B,X,O} = StridedBitPointer{N,C,B,R,X,NTuple{N,Int}}(p, dynamic_offsets(si))
@inline stridedpointer(p::Ptr{Bit}, si::StrideIndex{N,R,C,X,O}) where {N,R,C,X,O} = StridedBitPointer{N,C,0,R,X,NTuple{N,Int}}(p, dynamic_offsets(si))

@inline Base.pointer(p::Union{StridedPointer,StridedBitPointer}) = getfield(p, :p)
@inline ArrayInterface.StrideIndex(sptr::Union{StridedPointer,StridedBitPointer}) = getfield(sptr,:si)
@inline bytestrideindex(sptr::AbstractStridedPointer) = StrideIndex(sptr)

@inline bytestrides(si::StrideIndex) = _map(Base.Fix2(*,StaticInt{8}()), getfield(si,:strides))
@inline bytestrides(ptr::AbstractStridedPointer) = getfield(StrideIndex(ptr), :strides)
@inline Base.strides(ptr::AbstractStridedPointer) = getfield(StrideIndex(ptr), :strides)
@inline ArrayInterface.strides(ptr::AbstractStridedPointer) = getfield(StrideIndex(ptr), :strides)
@inline ArrayInterface.offsets(ptr::AbstractStridedPointer) = offsets(StrideIndex(ptr))
@inline ArrayInterface.contiguous_axis_indicator(ptr::AbstractStridedPointer{T,N,C}) where {T,N,C} = contiguous_axis_indicator(StaticInt{C}(), Val{N}())

@inline function similar_with_offset(sptr::AbstractStridedPointer{T,N,C,B,R}, ptr::Ptr, offset::Tuple) where {T,N,C,B,R}
  si = StrideIndex{N,R,C}(strides(sptr), offset)
  stridedpointer(ptr, si, contiguous_batch_size(sptr))
end
@inline function similar_no_offset(sptr::AbstractStridedPointer{T,N,C,B,R}, ptr::Ptr) where {T,N,C,B,R}
  si = StrideIndex{N,R,C}(strides(sptr), zerotuple(Val(N)))
  stridedpointer(ptr, si, contiguous_batch_size(sptr))
end

# There is probably a smarter way to do indexing adjustment here.
# The reasoning for the current approach of geping for Zero() on extracted inds
# and for offsets otherwise is best demonstrated witht his motivational example:
#
# A = OffsetArray(rand(10,11), 6:15, 5:15);
# for i in 6:15
    # s += A[i,i]
# end
# first access is at zero-based index
# (first(6:16) - ArrayInterface.offsets(a)[1]) * ArrayInterface.strides(A)[1] + (first(6:16) - ArrayInterface.offsets(a)[2]) * ArrayInterface.strides(A)[2]
# equal to
#  (6 - 6)*1 + (6 - 5)*10 = 10
# i.e., the 1-based index 11.
# So now we want to adjust the offsets and pointer's value
# so that indexing it with a single `i` (after summing strides) is correct.
# Moving to 0-based indexing is easiest for combining strides. So we gep by 0 on these inds.
# E.g., gep(stridedpointer(A), (0,0))
# ptr += (0 - 6)*1 + (0 - 5)*10 = -56
# new stride = 1 + 10 = 11
# new_offset = 0
# now if we access the 6th element
# (6 - new_offse) * new_stride
# ptr += (6 - 0) * 11 = 66
# cumulative:
# ptr = pointer(A) + 66 - 56  = pointer(A) + 10
# so initial load is of pointer(A) + 10 -> the 11th element w/ 1-based indexing

@inline stridedpointer(ptr::AbstractStridedPointer) = ptr
@inline stridedpointer_preserve(ptr::AbstractStridedPointer) = (ptr,nothing)

struct FastRange{T,F,S,O}# <: AbstractStridedPointer{T,1,1,0,(1,),Tuple{S},Tuple{O}}# <: AbstractRange{T}
  f::F
  s::S
  o::O
end
FastRange{T}(f::F,s::S) where {T<:Integer,F,S} = FastRange{T,Zero,S,F}(Zero(),s,f)
FastRange{T}(f::F,s::S,o::O) where {T,F,S,O} = FastRange{T,F,S,O}(f,s,o)

FastRange{T}(f::F,s::S) where {T<:FloatingTypes,F,S} = FastRange{T,F,S,Int}(f,s,0)
# FastRange{T}(f,s) where {T<:FloatingTypes} = FastRange{T}(f,s,fast_int64_to_double())
# FastRange{T}(f::F,s::S,::True) where {T<:FloatingTypes,F,S} = FastRange{T,F,S,Int}(f,s,0)
# FastRange{T}(f::F,s::S,::False) where {T<:FloatingTypes,F,S} = FastRange{T,F,S,Int32}(f,s,zero(Int32))

@inline function memory_reference(r::AbstractRange{T}) where {T}
  s = ArrayInterface.static_step(r)
  FastRange{T}(ArrayInterface.static_first(r) - s, s), nothing
end
@inline memory_reference(r::FastRange) = (r,nothing)
@inline bytestrides(::FastRange{T}) where {T} = (StaticInt(sizeof(T)),)
@inline ArrayInterface.offsets(::FastRange) = (One(),)
@inline val_stride_rank(::FastRange) = Val{(1,)}()
@inline val_dense_dims(::FastRange) = Val{(true,)}()
@inline ArrayInterface.contiguous_axis(::FastRange) = One()
@inline ArrayInterface.contiguous_batch_size(::FastRange) = Zero()

@inline stridedpointer(fr::FastRange, ::StrideIndex, ::StaticInt{0}) = fr
struct NoStrides end
@inline bytestrideindex(::FastRange) = NoStrides()
@inline ArrayInterface.offsets(::NoStrides) = NoStrides()
@inline reconstruct_ptr(r::FastRange{T}, o) where {T} = FastRange{T}(getfield(r,:f), getfield(r, :s), o)

@inline function zero_offsets(fr::FastRange{T,Zero}) where {T<:Integer}
  s = getfield(fr,:s)
  FastRange{T}(Zero(), s, getfield(fr,:o) + s)
end
@inline function zero_offsets(fr::FastRange{T}) where {T<:FloatingTypes}
  FastRange{T}(getfield(fr,:f), getfield(fr,:s), getfield(fr,:o) + 1)
end

# `FastRange{<:FloatingTypes}` must use an integer offset because `ptrforcomparison` needs to be exact/integral.

# discard unnueeded align/reg size info
@inline Base.eltype(::FastRange{T}) where {T} = T

# @inline Base.pointer(p::AbstractStridedPointer, i::Tuple) = pointer(p) + StrideIndex(p)[i...]
function reinterpret_strideindex_quote(
  sz_new::Int, sz_old::Int, C::Int, N::Int, curly::Expr
)
  if sz_old == sz_new
    # size_expr = :size_A
    bx_expr = :bx
  else
    @assert 1 ≤ C ≤ N
    # size_expr = Expr(:tuple)
    bx_expr = Expr(:tuple)
    for n ∈ 1:N
      # sz_n = Expr(:call, GlobalRef(Core,:getfield), :size_A, n, false)
      bx_n = Expr(:call, GlobalRef(Core,:getfield), :bx, n, false)
      if n ≠ C
        # push!(size_expr.args, sz_n)
        push!(bx_expr.args, bx_n)
      elseif sz_old > sz_new
        si = :(StaticInt{$(sz_old ÷ sz_new)}())
        # push!(size_expr.args, Expr(:call, :*, sz_n, si))
        push!(bx_expr.args, Expr(:call, :÷, bx_n, si))
      else
        si = :(StaticInt{$(sz_new ÷ sz_old)}())
        # push!(size_expr.args, Expr(:call, :÷, sz_n, si))
        push!(bx_expr.args, Expr(:call, :*, bx_n, si))
      end
    end
  end
  quote
    $(Expr(:meta,:inline))
    bx = si.strides
    $curly($bx_expr, si.offsets)
  end
end
@generated function reinterpret_strideindex(::Type{Tnew}, ::Type{Told}, si::StrideIndex{N,R,C}) where {Tnew,Told,N,R,C}
  # 1+2
  reinterpret_strideindex_quote(offsetsize(Tnew), offsetsize(Told), C, N, :(StrideIndex{$N,$R,$C}))
end
@generated function Base.reinterpret(::Type{Tnew}, sp::AbstractStridedPointer{Told,N,C,B,R}) where {Tnew,Told,N,C,B,R}
  # 1+2
  sz_old = offsetsize(Told)
  sz_new = offsetsize(Tnew)
  Bnew = ((B > 0) & (sz_new ≠ sz_old)) ? ((B * sz_old) ÷ sz_new) : B
  quote
    $(Expr(:meta,:inline))
    si = reinterpret_strideindex(Tnew, Told, StrideIndex(sp))
    stridedpointer(reinterpret(Ptr{$Tnew}, pointer(sp)), si, StaticInt{$Bnew}())
  end
end

@inline vload(::StaticInt{N}, args...) where {N} = StaticInt{N}()
@inline stridedpointer(::StaticInt{N}) where {N} = StaticInt{N}()
@inline zero_offsets(::StaticInt{N}) where {N} = StaticInt{N}()


text:
abstract type AbstractTCARule <: AbstractDCARule end

struct TCA{B,R,T} <: AbstractTCARule
    code::B
    codeset::R
    states::Int
    radius::T
end

"""
    TCA(code; states=2, radius=1)

Returns a ```TCA``` object given a specific code, number of states and radius. The ruleset for the rule is computed and
stored in the struct as well.
"""
function TCA(code;
    states=2,
    radius=1)

    codeset = tca_conversion(code, states, radius)
    TCA(code, codeset, states, radius)
end

function (tca::TCA)(starting_array)

    nextgen = tca_evolution(starting_array, tca.codeset, tca.states, tca.radius)
end

function tca_conversion(code, states, radius::Number)

    code_len = (2*radius+1)*states-2
    code_bin = parse.(Int, split(string(code, base=states), ""))
    code_bin = vcat(zeros(typeof(code_bin[1]), code_len-length(code_bin)), code_bin)
    reverse!(code_bin)
end

function tca_conversion(code, states, radius::Tuple)

    code_len = (sum(radius)+1)*states-2
    code_bin = parse.(Int, split(string(code, base=states), ""))
    code_bin = vcat(zeros(typeof(code_bin[1]), code_len-length(code_bin)), code_bin)
    reverse!(code_bin)
end

function tca_state_reader(neighborhood, codeset_len)

    mod1(sum(neighborhood)+1,codeset_len)
end

function tca_evolution(cell, codeset, states, radius::Number)

    neighborhood_size = radius*2+1
    output = zeros(length(cell))
    cell = vcat(cell[end-neighborhood_size÷2+1:end], cell, cell[1:neighborhood_size÷2])

    for i=1:length(cell)-neighborhood_size+1
        output[i] = codeset[tca_state_reader(cell[i:i+neighborhood_size-1], length(codeset))]
    end

    output
end

function tca_evolution(cell, codeset, states, radius::Tuple)

    neighborhood_size = sum(radius)+1
    output = zeros(length(cell))
    cell = vcat(cell[end-radius[1]+1:end], cell, cell[1:radius[2]])

    for i=1:length(cell)-neighborhood_size+1
        output[i] = codeset[tca_state_reader(cell[i:i+neighborhood_size-1], length(codeset))]
    end

    output
end


text:
using JuLIP

println("Test Direct use of ASECalculator")
at = bulk("Cu", cubic=true) * 2    # generate periodic Cu supercell
deleteat!(at, 1)                       # vacancy defect
try
   @pyimport ase.calculators.emt as emt   # import the EMT model
   calc = ASECalculator(emt.EMT())        # wrap it into a Julia Object
   @show energy(calc, at)                 # compute the energy
   @show maximum(norm.(forces(calc, at)))
   @test true
catch
   println("failed ASECalculator Test")
   @test false
end

# println("Test ASE's EMT calculator implementation (ASE.EMTCalculator)")
# emt = ASE.EMTCalculator()
# at = rattle!( set_pbc!( bulk("Cu", cubic=true) * 2, (true,false,false) ), 0.1 )
# @test JuLIP.Testing.fdtest(emt, at, verbose=true)
#
# println("Test ASE+JuLIP implementation of EMT")
# emtj = ASE.EMT.EMTCalculator(at)
#
# println("--------------------------------------------------")
# println(" EMT Consistency test: ")
# println("--------------------------------------------------")
# println(" E_ase - E_jl = ", energy(emt, at) - energy(emtj, at))
# println(" |Frc_ase - Frc_jl = ", maximum(norm.(forces(emt, at) - forces(emtj, at))))
# println("--------------------------------------------------")
# @test abs(energy(emt, at) - energy(emtj, at)) < 1e-10
# @test JuLIP.Testing.fdtest(emtj, at, verbose=true)


text:
"""
    DataPoint

Abstract type for KMeans poionts

Expect the following fields
# Fields
- `dimension`:
- `original`:
- `derived`:
"""
abstract type DataPoint end

struct SimpleDataPoint <: DataPoint
    dimension::Int64
    original::Vector{Float64}
    derived::Vector{Float64}

    function SimpleDataPoint(initial::Vector{Float64})
        dimension = length(initial)
        original = initial # reference to input array
        derived = copy(initial)

        new(dimension, original, derived)
    end
end

struct Governor <: DataPoint
    dimension::Int64
    original::Vector{Float64}
    derived::Vector{Float64}

    longitude::Float64
    age::Float64
    state::String

    function Governor(longitude, age, state)
        dimension = 2
        original = [longitude, age]
        derived = [longitude, age]
        new(dimension, original, derived, longitude, age, state)
    end

    function Governor(original::Vector{Float64})
        if 2 != length(original)
            throw(DomainError("Length of vector must be 2"))
        end
        dimension = 2
        derived = copy(original)
        new(2, original, derived, 0.0, 0.0, "")
    end
end

struct Album <: DataPoint
    dimension::Int64
    original::Vector{Float64}
    derived::Vector{Float64}

    name::String
    year::Int64
    length::Float64
    tracks::Float64

    function Album(name::String, year::Int64, length::Float64, tracks::Float64)
        dimension = 2
        original = [length, tracks]
        derived = [length, tracks]
        new(dimension, original, derived, name, year, length, tracks)
    end

    function Album(original::Vector{Float64})
        if 2 != length(original)
            throw(DomainError("Length of vector must be 2"))
        end
        dimension = 2
        derived = copy(original)
        new(2, original, derived, "", 0, 0.0, 0.0)
    end

end
function slice_name(pointvec::Vector{Album})
    vec = Vector{String}(undef, length(pointvec))

    for (i,p) in enumerate(pointvec)
        vec[i] = p.name
    end

    return vec
end

"""
    Base.:(==)(x::DataPoint, y::DataPoint)

Check `derived` vector for equality
"""
function Base.:(==)(x::DataPoint, y::DataPoint)
    return x.derived == y.derived
end

"""
    set_original!(x::P, initial) where {P<:DataPoint}

Convenient function to set `original` and `derived` to the same
values.
"""
function set_original!(x::P, initial) where {P<:DataPoint}
    if x.dimension != length(initial)
        throws(DomainError("Dimension mismatching"))
    end

    x.original .= initial
    x.derived .= initial
end

function distance(x::P, y::P) where {P<:DataPoint}
    if x.dimension != y.dimension
        throw(DomainError("Mis-matching dimensions: $(x.dimension), $(y.dimension)"))
    end

    diff2 = @. (x.derived - y.derived) ^ 2

    return sqrt(sum(diff2))
end

"""
    Cluster

A cluster of data points

# Fields
- `pointvec`: points
- `centroid`: centroid of points
"""
struct Cluster{P<:DataPoint}
    pointvec::Vector{P}
    centroid::P
end

"""
    KMeans{P<:DataPoint}

K-Mean data structure with points and clusters

# Fields
- `pointvec`: Vector of DataPoint
- `clustervec`: Vector of Cluster

"""
struct KMeans{P<:DataPoint}
    pointvec::Vector{P}
    clustervec::Vector{Cluster{P}}
end

"""
    KMeans(k::Int64, pointvec::Vector{P}) where {P<:DataPoint}

Creates a `KMeans` with `k` clusters and random centroids using the
input points.
"""
function KMeans(k::Int64, pointvec::Vector{P}) where {P<:DataPoint}
    if k < 1
        throw(DomainError("k must be >= 1"))
    end

    clustervec = Vector{Cluster{P}}()
    me = KMeans(pointvec, clustervec)

    # initialize with random clusters
    #normalize_zscore!(me.pointvec)
    zscore_derived!(me.pointvec)

    #dimension_count = pointvec[1].dimension
    for i in 1:k
        rand_point = random_point(pointvec)
        cluster = Cluster(Vector{P}(), rand_point)
        push!(me.clustervec,cluster)
    end

    return me
end

"""
    slice_derived(pointvec::Vector{P}, dimension) where {P<:DataPoint}

Cut a slice out of `DataPoint` at the specified dimension.

`dimensionSlice` in Java
"""
function slice_derived(pointvec::Vector{P}, dimension) where {P<:DataPoint}
    derivedvec = Vector{Float64}(undef, length(pointvec))

    for (i,p) in enumerate(pointvec)
        derivedvec[i] = p.derived[dimension]
    end

    return derivedvec
end

"""
    slice_original(pointvec::Vector{P}, dimension) where {P<:DataPoint}

Get `original[dimension]` from each point and put into a vector.

Similar to `dimensionSlice` in Java.
"""
function slice_original(pointvec::Vector{P}, dimension) where {P<:DataPoint}
    originalvec = Vector{Float64}(undef, length(pointvec))

    for (i,p) in enumerate(pointvec)
        originalvec[i] = p.original[dimension]
    end

    return originalvec
end

function slice_state(pointvec::Vector{Governor})
    statevec = Vector{String}(undef, length(pointvec))

    for (i,p) in enumerate(pointvec)
        statevec[i] = p.state
    end

    return statevec
end

"""
    normalize_zscore!(pointvec::Vector{P}) where {P<:DataPoint}

Assign zscore of `derived` to `derived`

Relys on `P` has just been initialized so `derived` equals to `original`.
Do not think this is a good approch.  See `zscore_derived()`.

Create `derived` data that is normalized using zscore.
The is the `zScoreNormalize` in Java.
"""
function normalize_zscore!(pointvec::Vector{P}) where {P<:DataPoint}
    zscored = Vector{Vector{Float64}}()
    for i in 1:length(pointvec)
        push!(zscored, Vector{Float64}())
    end

    dimension = pointvec[1].dimension
    for dim_ind = 1:dimension
        dimension_slice = slice_derived(pointvec, dimension)
        mu = mean(dimension_slice)
        sigma = std(dimension_slice, corrected=false)
        zscores = zscore(dimension_slice, mu, sigma)
        for ind in 1:length(zscores)
            push!(zscored[ind], zscores[ind])
        end

    end

    for i in 1:length(pointvec)
        pointvec[i].derived .= zscored[i]
    end
end

"""
    zscore_derived!(pointvec::Vector{P}) where {P<:DataPoint}

Assign zscore values of `original` to `derived`

This resembles the `zScoreNormalize` in Java but not exactly the same.
"""
function zscore_derived!(pointvec::Vector{P}) where {P<:DataPoint}
    zscored = Vector{Vector{Float64}}()
    for i in 1:length(pointvec)
        push!(zscored, Vector{Float64}())
    end

    dimension = pointvec[1].dimension
    for dim_ind = 1:dimension
        dimension_slice = slice_original(pointvec, dim_ind)
        mu = mean(dimension_slice)
        sigma = std(dimension_slice, corrected=false)
        zscores = zscore(dimension_slice, mu, sigma)
        for ind in 1:length(zscores)
            push!(zscored[ind], zscores[ind])
        end

    end

    for i in 1:length(pointvec)
        pointvec[i].derived .= zscored[i]
    end
end

"""
    random_point(pointvec::Vector{P}) where {P<:DataPoint}

Create a random point using the `derived` data.  The point
is within the min/max of the derived data.
"""
function random_point(pointvec::Vector{P}) where {P<:DataPoint}
    rand_derived = Vector{Float64}()
    dimension_count = pointvec[1].dimension
    for i in 1:dimension_count
        values = slice_derived(pointvec, i)
        Max = maximum(values)
        Min = minimum(values)
        rand_value = Min + rand() * (Max-Min)
        push!(rand_derived, rand_value)
    end

    return P(rand_derived)
end

"""
    assign_clusters!(kmeans::KMeans)

Update the points of the cluster of `KMeans` based on
the centroid of the cluster of `KMeans`.  That is update
`kmeans.clustervec.pointvec` based on `kmeans.clustervec.centroid`
"""
function assign_clusters!(kmeans::KMeans)
    for cluster in kmeans.clustervec
        resize!(cluster.pointvec,0) # clear
    end

    for point in kmeans.pointvec
        lowest_distance = typemax(Float64)
        cloest_cluster = kmeans.clustervec[1]
        for cluster in kmeans.clustervec
            centroid_distance = distance(point, cluster.centroid)
            if centroid_distance < lowest_distance
                lowest_distance = centroid_distance
                cloest_cluster = cluster
            end
        end

        # pointvec was cleared outside
        push!(cloest_cluster.pointvec, point)
    end
end

"""
    centroids(kmeans::KMeans)

Get current centroids of `KMeans`

This can be used to check if successive centroids changed to
determine convergence.
"""
function centroids(kmeans::KMeans{P})::Vector{P} where {P<:DataPoint}
    [c.centroid for c in kmeans.clustervec]
end

"""
    generate_centroids!(kmeans::KMeans)

Update the centroid of each cluster based on points of cluster.
"""
function generate_centroids!(kmeans::KMeans)
    for cluster in kmeans.clustervec
        if isempty(cluster.pointvec)
            continue
        end

        means = Vector{Float64}()
        dimension_count = cluster.pointvec[1].dimension
        for dimension in 1:dimension_count

            derived = slice_derived(cluster.pointvec, dimension)
            derived_mean = mean(derived)
            push!(means, derived_mean)
        end

        set_original!(cluster.centroid, means)
    end
end

"""
    run!(kmeans::KMeans, maxIt::Int64)

Run k-mean clustering

Return number of iterations
"""
function run!(kmeans::KMeans, maxIt::Int64)
    for it in 1:maxIt
        # Moved to assign_clusters!
        # for cluster in kmeans.clustervec
        #     resize!(cluster.pointvec,0) # clear
        # end

        assign_clusters!(kmeans)
        old_centroids = centroids(kmeans)
        generate_centroids!(kmeans)
        if old_centroids == centroids(kmeans)
            #println("Converged after $it iterations.")
            return it
        end
    end
end


text:
__precompile__

module ColloidFluctuatingField
export time_propagation, para
using FFTW
using Random
using Parameters

include("helper.jl")
include("modelA_fourier_galerkin.jl")
include("NewtonImp.jl")
include("solver.jl")
include("driven_cond.jl")

@with_kw struct para
    Nt::Int64
    Δt::Float64
    L::Float64
    N::Int64
    Φ::Array
    x::Array{Float64,1}
    ẋ::Array{Float64,1}
    para_Φ::Array
    para_x::Array
    para_V::Array
    times_save::Array{Int64,1}=reverse([i for i = 1:Nt])
end

function check_method(method::String, dic::Dict, st::String)
    if haskey(dic,method)
        ans = dic[method]
        if ans === nothing
            error("Please provide costum method for time propagation of" *st)
        end
    else
        error("Method not found for "*st)
    end
    return ans
end

function time_propagation(p::para, method_Φ::String,
    method_x::String, method_V::String, custom_Φ = nothing, custom_x = nothing,
    custom_V = nothing)
    Φ_methods = Dict("ModelA-Galerkin" =>Φ_galerk_modelA,
                    "Driven-Condensate" => Φ_galerk_driven_cond, "custom" => custom_Φ)
    x_methods = Dict("Newton"=> NewtonImp, "custom" =>custom_x)
    V_methods = Dict("Gaussian" => V_gaussian, "custom" =>custom_V)
    ɸ_fct = check_method(method_Φ, Φ_methods,"Φ")
    V_fct = check_method(method_V, V_methods,"V")
    x_fct = check_method(method_x, x_methods, "x")
    @unpack Nt, Δt, L, N, Φ, x, ẋ, para_Φ, para_x, para_V, times_save = p
    ɸ_r, x_r, ẋ_r, ẍ_r = solver(Nt, Δt, L, N, Φ, x, ẋ, para_Φ, para_x, para_V, ɸ_fct, x_fct, V_fct, times_save)
    return ɸ_r, x_r, ẋ_r, ẍ_r, p
end

end

max_stars_repo_name: phyjonas/ColloidFluctuatingField


text:
import Base:show, size, summary

const show_os = 12
const unindexed_fields = (:v, :n)

summary(B::BeamData) = string(typeof(B), " with ", length(B.t), " beampower",
                             size(B.bp,2) == 1 ? "" : "s")

function show(io::IO, B::BeamData)
  W = max(80, displaysize(io)[2]) - show_os
  w = min(W, 35)
  nc = size(B.bp,1) == 0 ? 0 : 1
  N = min(nc, div(W-1, w))
  M = min(N+1, nc)
  bplen,numbp = size(B.bp)
  println(io, "BeamData with ", length(B.t), " beampowers")
  Field = fieldnames(BeamData)
  for f in Field
    if (f in unindexed_fields) == false
      targ = getfield(B, f)
      t = typeof(targ)
      fstr = uppercase(String(f))
      print(io, lpad(fstr, show_os-2), ": ")
      if f == :notes || f == :misc
        SeisIO.show_str(io, String[string(length(targ), " entries")], w, N)
      elseif (t <: AbstractFloat || t <: InstrumentPosition || t<: InstrumentResponse)
        println(io, repr("text/plain", targ, context=:compact=>true))
      elseif f == :t
        if length(B.t) > 0
          SeisIO.show_str(io, String[timestamp(B.t[1]), " (", string(length(B.t)), " BeamPowers)"], w, N)
        else
          SeisIO.show_str(io, String[repr("text/plain", targ, context=:compact=>true)], w, N)
        end
      elseif f == :bp
        SeisIO.show_str(io, String[summary(targ)], w, N)
      else
        SeisIO.show_str(io, String[repr("text/plain", targ, context=:compact=>true)], w, N)
      end
    end
  end
  return nothing
end
show(B::BeamData) = show(stdout, B)

max_stars_repo_name: kura-okubo/SeisBeamforming

max_issues_repo_name: kura-okubo/SeisBeamforming

max_forks_repo_name: kura-okubo/SeisBeamforming


text:
# Plots of quadrature geometry

push!(LOAD_PATH, string(pwd(),"/src"))

using Revise
using PyPlot
import IterativeSolvers
import LinearMaps

import ModifiedStokesSolver
import AnalyticDomains

mss = ModifiedStokesSolver

# Params
alpha = 200

# Discretize
numpanels = 20
panelorder = 16
curve = AnalyticDomains.starfish(n_arms = 3, amplitude=0.1)
dcurve = mss.discretize(curve, numpanels, panelorder)

N = dcurve.numpoints
panel_length = sum(dcurve.dS[1:panelorder])

n_sub = 32
dt_max = mss.LARGEALPHA_LIMIT_16*2/(alpha*panel_length)
glpoints1, glweights1 = mss.gausslegendre(panelorder)
glpoints2, glweights2 = mss.gausslegendre(n_sub)

######################## Near boundary point
panel_idx = 1
point_idx = 6

t1 = dcurve.t_edges[:,panel_idx]
idx = panelorder*(panel_idx-1) + (1:panelorder)
xpanel = dcurve.points[1,idx]
ypanel = dcurve.points[2,idx]
zj = xpanel + 1im*ypanel
nj = dcurve.normals[1,idx] + 1im*dcurve.normals[2,idx]

h1 = t1[2]-t1[1]
@show t = glpoints1[point_idx] + 0.08im;
ztmp = curve.tau( t1[1] + (t+1)*h1/2 )
zt = [real(ztmp), imag(ztmp)]

coeffs = CurveDiscretization.map_panels(dcurve)
zcpx = zt[1] + 1im*zt[2]
troot, _, converged =
    CurveDiscretization.invert_map(dcurve, coeffs, panel_idx, zcpx,
                                   fall_back_to_initial=false)
@show troot

paperplots = true
if paperplots
    close("all")
    figure(1, figsize=(3,3))
else
    figure(1)
    clf()
end

plot(dcurve.points[1,:], dcurve.points[2,:], "-k", label=L"\Gamma")

cpx2vec(z) = [real(z),imag(z)]
function plot_edge(x, n, h, style="b")
    plot(x[1]+n[1]*[-1,1]*h, x[2]+n[2]*[-1,1]*h, style)
end

h1 = 0.06*panel_length
h2 = 0.03*panel_length

for i=1:numpanels
    x = dcurve.edges[1:2,i]
    n = dcurve.n_edges[1:2,i]
    plot_edge(x, n, h1)
end

plot(zt[1], zt[2], ".", label=L"x \in \Omega")

# Plot subintervals
bcweights = mss.bclag_interp_weights(glpoints1)
@show intervals = mss.subdivide_interval_with_bisection(dt_max, troot, n_sub)
#intervals = mss.subdivide_interval(dt_max, troot)
for idx = 1:length(intervals)-1
    ta = intervals[idx]
    tb = intervals[idx+1]
    troot_sec = (troot-ta)*2/(tb-ta) - 1.0
    # New source points, in original frame
    tsrc = ta + (1+glpoints2)*(tb-ta)/2
    # Interpolation to new source points
    Psec = mss.bclag_interp_matrix(glpoints1, tsrc, bcweights)
    zj_sec = Psec*zj
    za_sec = mss.bclag_interp_eval(glpoints1, zj, ta, bcweights)
    zb_sec = mss.bclag_interp_eval(glpoints1, zj, tb, bcweights)
    nj_sec = Psec*nj
    na_sec = mss.bclag_interp_eval(glpoints1, nj, ta, bcweights)
    nb_sec = mss.bclag_interp_eval(glpoints1, nj, tb, bcweights)

    if idx>1
        plot_edge(cpx2vec(za_sec), cpx2vec(na_sec), h2, "r")
    end
    if idx<length(intervals)-1
        plot_edge(cpx2vec(zb_sec), cpx2vec(nb_sec), h2, "r")
    end
end

xticks([])
yticks([])
minorticks_off()
axis("equal")
axis((0.9153349805587189, 1.2063410174030949, -0.02194527272134522, 0.3516612032549363))
legend()

function write_fig(num, name)
    figure(num)
    savefig(name)
    println("Saved $name")
end

filename = joinpath(dirname(@__DIR__),"figs", "quad_geo_near.pdf")
write_fig(1, filename)
run(`cp $filename ../docs/paper/figs`)

######################## On boundary point
zt = dcurve.points[:,point_idx]

coeffs = CurveDiscretization.map_panels(dcurve)
zcpx = zt[1] + 1im*zt[2]
troot = glpoints1[point_idx]

paperplots = true
if paperplots
    figure(2, figsize=(3,3))
else
    figure(2)
    clf()
end

plot(dcurve.points[1,:], dcurve.points[2,:], "-k", label=L"\Gamma")

for i=1:numpanels
    x = dcurve.edges[1:2,i]
    n = dcurve.n_edges[1:2,i]
    plot_edge(x, n, h1)
end

plot(zt[1], zt[2], ".", label=L"x \in \Gamma")

# Plot subintervals
intervals = mss.subdivide_interval_with_bisection(dt_max, troot, 32)
#intervals = mss.subdivide_interval(dt_max, troot)
for idx = 2:length(intervals)-2
    ta = intervals[idx]
    tb = intervals[idx+1]
    troot_sec = (troot-ta)*2/(tb-ta) - 1.0
    # New source points, in original frame
    tsrc = ta + (1+glpoints2)*(tb-ta)/2
    # Interpolation to new source points
    Psec = mss.bclag_interp_matrix(glpoints1, tsrc, bcweights)
    zj_sec = Psec*zj
    za_sec = mss.bclag_interp_eval(glpoints1, zj, ta, bcweights)
    zb_sec = mss.bclag_interp_eval(glpoints1, zj, tb, bcweights)
    nj_sec = Psec*nj
    na_sec = mss.bclag_interp_eval(glpoints1, nj, ta, bcweights)
    nb_sec = mss.bclag_interp_eval(glpoints1, nj, tb, bcweights)

    plot_edge(cpx2vec(za_sec), cpx2vec(na_sec), h2, "r")
    plot_edge(cpx2vec(zb_sec), cpx2vec(nb_sec), h2, "r")
end

xticks([])
yticks([])
minorticks_off()
axis("equal")
axis((0.9153349805587189, 1.2063410174030949, -0.02194527272134522, 0.3516612032549363))
legend()

filename = joinpath(dirname(@__DIR__),"figs", "quad_geo_on.pdf")
write_fig(2, filename)
run(`cp $filename ../docs/paper/figs`)


text:
using PyCall
using LinearAlgebra

pushfirst!(PyVector(pyimport("sys")."path"), "")
pymatrix = pyimport("matrix")

matArray = pymatrix.read_csv("matrix.csv")
println("Input Matrix:")
println(matArray)
println("Determinant:")
println(det(matArray))
println("Trace:")
println(tr(matArray))
println("Square:")
squaredArray = matArray * matArray
println(squaredArray)
destroyedArray = squaredArray * squaredArray * det(squaredArray)
println("Assortment of various matrix operations:\n", destroyedArray)
println("Wow! That was fast!")
println("Writing to file...")
pymatrix.write_csv(destroyedArray, "output.csv")

max_stars_repo_name: calvang/pyjulia-experiments

max_issues_repo_name: calvang/pyjulia-experiments

max_forks_repo_name: calvang/pyjulia-experiments


text:
    using FileIO
using LinearAlgebra
using GLMakie, AbstractPlotting
AbstractPlotting.__init__()

dt = 0.0001
N_Schritt = 50

dt = 60*60. / N_Schritt

AU = 149.6e9 # m
Masse_Erde = 5.972e24 # kg
Radius_Erde = 6378_000f0 # m

mutable struct Körper
    Position
    Geschinwigkeit
    Radius
    Masse
    Bild
    Object
    Attribute
end

function Körper(scene,   Position,   Geschinwigkeit,   Radius,  Masse,  Bild)
    Object = mesh!(scene,Sphere(Point3f0(0), Float32(Radius)), color = load(Bild), shading = false)[end]
    translate!(Object, Position);
    return Körper(Position,   Geschinwigkeit,   Radius,  Masse, Bild, Object, Dict())
end

position(K) = translation(K.Object)[]

scene = Scene(show_axis = false, center = false)

f_sun = 20*Radius_Erde
f = 1000 * Radius_Erde
#f = 200 * Radius_Erde

System = [Körper(scene,Point3f0(0,0,0), Point3f0(0,0,0), 109*f_sun, 1.989e30, "Sonne.jpg"),
          Körper(scene,Point3f0(0.387 * AU,0,0), Point3f0(0,47.4e3,0), 0.382 * f, 0.055 * Masse_Erde, "Merkur.png"),
          Körper(scene,Point3f0(0.723 * AU,0,0), Point3f0(0,35e3,0), 0.949 * f, 0.8150 * Masse_Erde, "Venus.jpg"),
          Körper(scene,Point3f0(AU,0,0), Point3f0(0,29.780e3,0), f, Masse_Erde, "Erde.png"),
          Körper(scene,Point3f0(1.524  * AU,0,0), Point3f0(0,24e3,0), 0.532 * f, 0.64171e24, "Mars.jpg"),
          Körper(scene,Point3f0(5.2 * AU,0,0), Point3f0(0,13.1e3,0), 11.19 * f, 1.898e27, "Jupiter.jpg"),
          Körper(scene,Point3f0(AU+405400e3,0,0), Point3f0(0,1.022e3 + 29.780e3 ,0), 0.2727 * f, 0.012 * Masse_Erde, "Mond.jpg")
          ]

#update_cam!(scene, Vec3f0(-10e10, 10e10, 10e10), Vec3f0(0))
update_cam!(scene, Vec3f0(0e10, 0e10, 10e10), Vec3f0(0), Vec3f0(0, 1, 0) )
scene.center = false # prevent to recenter on display

display(scene)

Zeit = Ref(0.)
Kräfte = zeros(length(System),3)
Positionen = zeros(length(System),3)
Geschwindigkeiten = zeros(length(System),3)
nPositionen = zeros(length(System),3)
nGeschwindigkeiten = zeros(length(System),3)

for i = 1:length(System)
    Geschwindigkeiten[i,:] = System[i].Geschinwigkeit;
end

function neuePosition(System,dt,Positionen,Geschwindigkeiten,Kräfte,nPositionen,nGeschwindigkeiten)
    G = 6.67430e-11 # m³/(kg s²)
    Kräfte .= 0

    for i = 1:length(System)
        for j = 1:length(System)
            if i ≠ j
                r = Positionen[i,:] - Positionen[j,:]
                Kraft = G * System[i].Masse * System[j].Masse * r ./ norm(r)^3
                Kräfte[i,:] = Kräfte[i,:] - Kraft
                Kräfte[j,:] = Kräfte[j,:] + Kraft
            end
        end
    end

    for i = 1:length(System)
        nGeschwindigkeiten[i,:] = Geschwindigkeiten[i,:] + Kräfte[i,:]/System[i].Masse * dt
        nPositionen[i,:] = Positionen[i,:] + Geschwindigkeiten[i,:] * dt
    end

end

function Zeitschritt(System,Zeit,Kräfte)

    for i = 1:length(System)
        Positionen[i,:] = translation(System[i].Object)[];
    end

    for n = 1:N_Schritt
        neuePosition(System,dt,Positionen,Geschwindigkeiten,Kräfte,nPositionen,nGeschwindigkeiten)
        @. nPositionen = (Positionen+nPositionen) * 0.5
        @. nGeschwindigkeiten = (Geschwindigkeiten+nGeschwindigkeiten) * 0.5
        neuePosition(System,dt,nPositionen,nGeschwindigkeiten,Kräfte,Positionen,Geschwindigkeiten)
    end

    for i = 1:length(System)
        t = translation(System[i].Object);
        t[] = Positionen[i,:]
    end

end

function anim(scene,System,Zeit,Kräfte)
    while isopen(scene)
        Zeit[] += 0.1

	    Zeitschritt(System,Zeit,Kräfte)
	    sleep(0.01)
    end
end

#task = @async anim(scene,System,Zeit,Kräfte)

anim(scene,System,Zeit,Kräfte)
stop() = schedule(task,InterruptException,error = true)

#=
close = false
while !close
global close
h = on(scene.events.window_open) do val
if !val
close = true
end
end
sleep(0.1)
end
=#

#Adrian


text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl

combined predict and center stepper
=#

mutable struct CombinedStepper{T <: Real} <: Stepper{T}
    shift_sched::Int
    searcher_options

    prev_alpha::T
    rhs::Point{T}
    dir::Point{T}
    temp::Point{T}
    dir_cent::Point{T}
    dir_centadj::Point{T}
    dir_pred::Point{T}
    dir_predadj::Point{T}
    dir_temp::Vector{T}

    searcher::StepSearcher{T}
    cent_only::Bool
    unadj_only::Bool

    function CombinedStepper{T}(;
        shift_sched::Int = 0,
        searcher_options...
        ) where {T <: Real}
        stepper = new{T}()
        stepper.shift_sched = shift_sched
        stepper.searcher_options = searcher_options
        return stepper
    end
end

function load(stepper::CombinedStepper{T}, solver::Solver{T}) where {T <: Real}
    model = solver.model

    stepper.prev_alpha = one(T)
    stepper.rhs = Point(model)
    stepper.dir = Point(model)
    stepper.temp = Point(model)
    stepper.dir_cent = Point(model, ztsk_only = true)
    stepper.dir_pred = Point(model, ztsk_only = true)
    stepper.dir_centadj = Point(model, ztsk_only = true)
    stepper.dir_predadj = Point(model, ztsk_only = true)
    stepper.dir_temp = zeros(T, length(stepper.rhs.vec))

    stepper.searcher = StepSearcher{T}(model; stepper.searcher_options...)
    stepper.unadj_only = stepper.cent_only = false

    return stepper
end

function step(stepper::CombinedStepper{T}, solver::Solver{T}) where {T <: Real}
    point = solver.point
    model = solver.model
    rhs = stepper.rhs
    dir = stepper.dir
    dir_cent = stepper.dir_cent
    dir_centadj = stepper.dir_centadj
    dir_pred = stepper.dir_pred
    dir_predadj = stepper.dir_predadj

    # update linear system solver factorization
    solver.time_upsys += @elapsed update_lhs(solver.syssolver, solver)

    # calculate centering direction and adjustment
    solver.time_uprhs += @elapsed update_rhs_cent(solver, rhs)
    solver.time_getdir += @elapsed get_directions(stepper, solver)
    copyto!(dir_cent.vec, dir.vec)
    solver.time_uprhs += @elapsed update_rhs_centadj(solver, rhs, dir)
    solver.time_getdir += @elapsed get_directions(stepper, solver)
    copyto!(dir_centadj.vec, dir.vec)

    # calculate affine/prediction direction and adjustment
    solver.time_uprhs += @elapsed update_rhs_pred(solver, rhs)
    solver.time_getdir += @elapsed get_directions(stepper, solver)
    copyto!(dir_pred.vec, dir.vec)
    solver.time_uprhs += @elapsed update_rhs_predadj(solver, rhs, dir)
    solver.time_getdir += @elapsed get_directions(stepper, solver)
    copyto!(dir_predadj.vec, dir.vec)

    # search with combined directions and adjustments
    stepper.unadj_only = stepper.cent_only = false
    solver.time_search += @elapsed alpha = search_alpha(point, model, stepper)

    if iszero(alpha)
        # recover
        solver.verbose && println("trying combined without adjustment")
        stepper.unadj_only = true
        solver.time_search += @elapsed alpha = search_alpha(point, model, stepper)

        if iszero(alpha)
            solver.verbose && println("trying centering with adjustment")
            stepper.cent_only = true
            stepper.unadj_only = false
            solver.time_search += @elapsed alpha =
                search_alpha(point, model, stepper)

            if iszero(alpha)
                solver.verbose && println("trying centering without adjustment")
                stepper.unadj_only = true
                solver.time_search += @elapsed alpha =
                    search_alpha(point, model, stepper)

                if iszero(alpha)
                    @warn("cannot step in centering direction")
                    solver.status = NumericalFailure
                    stepper.prev_alpha = alpha
                    return false
                end
            end
        end
    end

    # step
    update_stepper_points(alpha, point, stepper, false)
    stepper.prev_alpha = alpha

    return true
end

expect_improvement(stepper::CombinedStepper) = true

function update_stepper_points(
    alpha::T,
    point::Point{T},
    stepper::CombinedStepper{T},
    ztsk_only::Bool,
    ) where {T <: Real}
    if ztsk_only
        cand = stepper.temp.ztsk
        copyto!(cand, point.ztsk)
        dir_cent = stepper.dir_cent.ztsk
        dir_pred = stepper.dir_pred.ztsk
    else
        cand = point.vec
        dir_cent = stepper.dir_cent.vec
        dir_pred = stepper.dir_pred.vec
    end

    if stepper.unadj_only
        # no adjustment
        if stepper.cent_only
            # centering
            @. cand += alpha * dir_cent
        else
            # combined
            alpha_m1 = 1 - alpha
            @. cand += alpha * dir_pred + alpha_m1 * dir_cent
        end
    else
        # adjustment
        dir_centadj = (ztsk_only ? stepper.dir_centadj.ztsk :
            stepper.dir_centadj.vec)
        alpha_sqr = abs2(alpha)
        if stepper.cent_only
            # centering
            @. cand += alpha * dir_cent + alpha_sqr * dir_centadj
        else
            # combined
            dir_predadj = (ztsk_only ? stepper.dir_predadj.ztsk :
                stepper.dir_predadj.vec)
            alpha_m1 = 1 - alpha
            alpha_m1sqr = abs2(alpha_m1)
            @. cand += alpha * dir_pred + alpha_sqr * dir_predadj +
                alpha_m1 * dir_cent + alpha_m1sqr * dir_centadj
        end
    end

    return
end

function start_sched(stepper::CombinedStepper, searcher::StepSearcher)
    (stepper.shift_sched <= 0) && return 1
    return max(1, searcher.prev_sched - stepper.shift_sched)
end

print_header_more(stepper::CombinedStepper, solver::Solver) =
    @printf("%5s %9s", "step", "alpha")

function print_iteration_more(stepper::CombinedStepper, solver::Solver)
    if stepper.cent_only
        step = (stepper.unadj_only ? "cent" : "ce-a")
    else
        step = (stepper.unadj_only ? "comb" : "co-a")
    end
    @printf("%5s %9.2e", step, stepper.prev_alpha)
    return
end


text:
type Index
    lookup::Dict{Symbol, Int}
    names::Vector{Symbol}
end

function Index(names::Vector{Symbol})
    lookup = Dict{Symbol, Int}(zip(names, 1:length(names)))
    Index(lookup, names)
end

type GroupBy
    index::Index
    gvars::Vector{Int}          # indexes of vars to group by
    groups::Dict{Any, Vector{VarStats}}
    groupers::Dict{Symbol, Function}
end

function GroupBy(names::Vector{Symbol}, groupers::Dict{Symbol, Function})
    index = Index(names)
    gnames = filter(x -> in(x, keys(groupers)), names)
    gvars = [index.lookup[name] for name in gnames]
    return GroupBy(index, gvars, Dict{Any, Vector{VarStats}}(), groupers)
end

function GroupBy(names::Vector{Symbol}, groupby::Vector{Symbol})
    return GroupBy(names, Dict([(name, identity) for name in groupby]))
end

show(io::IO, grb::GroupBy) =
    print(io, "GroupBy($(grb.index.names),$(keys(grb.groupers)))")

function applygroupers{T<:Real}(grb::GroupBy, values::Vector{T})
    @assert length(grb.index.names) == length(values) "number of names != number of values"
    grouped = copy(values)
    for name in keys(grb.groupers)
        idx = grb.index.lookup[name]
        grouped[idx] = grb.groupers[name](grouped[idx])
    end
    return grouped
end

function update!{T<:Real}(grb::GroupBy, values::Vector{T})
    groupedvalues = applygroupers(grb, values)
    key = groupedvalues[grb.gvars]
    if !haskey(grb.groups, key)
        # init new varstats group to the length of first incoming values
        grb.groups[key] = [VarStats() for i=1:length(groupedvalues)]
    end
    updatestats!(grb.groups[key], groupedvalues)
end

function minimum(grb::GroupBy, name)
    idx = grb.index.lookup[name]
    minimum([minimum(grb.groups[k][idx]) for k in keys(grb.groups)])
end

function maximum(grb::GroupBy, name)
    idx = grb.index.lookup[name]
    maximum([maximum(grb.groups[k][idx]) for k in keys(grb.groups)])
end

function sum(grb::GroupBy, name::Symbol)
    idx = grb.index.lookup[name]
    sum([sum(grb.groups[k][idx]) for k in keys(grb.groups)])
end

function count(grb::GroupBy, name::Symbol)
    idx = grb.index.lookup[name]
    sum([count(grb.groups[k][idx]) for k in keys(grb.groups)])
end

function mean(grb::GroupBy, name::Symbol)
    sum(grb, name) / count(grb, name)
end

function hist_from_counts{T<:Real}(counts::Dict{T, Int})
    sum_counts = Dict(zip(keys(counts), zeros(Int, length(counts))))
    for (v, c) in counts
        sum_counts[v] += c
    end
    sorted_counts = sort([(v, c) for (v, c) in sum_counts], by=b->b[1])
    bins = [v for (v, c) in sorted_counts]
    counts = [c for (v, c) in sorted_counts]
    return bins, counts
end

function hist(grb::GroupBy, name::Symbol)
    idx = grb.index.lookup[name]
    @assert in(idx, grb.gvars) "Can only build histogram on grouper variables"
    gidx = findin(grb.gvars, idx)[1]
    counts = Dict{Float64, Int}([(k[gidx], count(vss[idx]))
                                 for (k, vss) in grb.groups])
    return hist_from_counts(counts)
end

function quantile(grb::GroupBy, name::Symbol, q::Real)
    @assert (0 <= q <= 1) "Quantile is $q, should be in [0..1]"
    (bins, h) = hist(grb, name)
    qs = q * sum(h)
    s = 0
    i = 0
    while s < qs
        i += 1
        s += h[i]
    end
    return if i > 0 bins[i] else 1. end
end

# groupers

@doc """
Create grouper that adjusts value to the nearest multiple of step.
Usage:
    grouper = stepgrouper(5)
    grouper(75) ==> 75
    grouper(78) ==> 80
    grouper(80) ==> 80
""" ->
function stepgrouper(step::Real)
    function grouper(value)
        round(value / step) * step
    end
    return grouper
end


text:
using Base: @deprecate

@deprecate Distinct() Unique()
@deprecate TeeZip(xf) ZipSource(xf) false

# It wasn't a public API but to play on the safe side:
@deprecate induction(itr) extract_transducer(itr) false


text:
include("opt_pid_data.jl")
include("opt_pid_work.jl")
module Opt_PID

using Opt_PID_data: My_Eval, create_My_Eval, set_rhs!,  condEntropy
using Mosek
using MathProgBase

struct PID_Results
    support  :: Array{Tuple{Int,Int,Int}}    # support of the following two vectors
    q        :: Vector{Float64}              # argmax H(X|YZ)
    # ∇        :: Vector{Float64}              # (super-)gradient of H(X|YZ) in p
    λ        :: Vector{Float64}              # (super-)gradient of H(X|YZ) in p
    p        :: Vector{Float64}              # current iterate (3-dim array!)
end

struct PID_Data
    bits     :: Float64
    myeval   :: My_Eval

    results  :: PID_Results

    # Storage --- this things are only here because I don't want to alloc mem for them in every iteration
    _space   :: Array{Float64,3} # 3D-storage for p
    _lb      :: Vector{Float64}
    _ub      :: Vector{Float64}
    _l       :: Vector{Float64}
    _u       :: Vector{Float64}
end

function init_pid_data(n_X::Int, n_Y::Int, n_Z::Int, their_support::Array{Tuple{Int,Int,Int}} ; tmpFloat::DataType=BigFloat) :: PID_Data
    pdummy = zeros(Float64,n_X,n_Y,n_Z)
    for xyz in their_support
        x,y,z = xyz
        @assert 1 ≤ x ≤ n_X "x out of range"
        @assert 1 ≤ y ≤ n_Y "y out of range"
        @assert 1 ≤ z ≤ n_Z "z out of range"

        pdummy[ x,y,z ] = 1.
    end

    myeval = create_My_Eval(pdummy,tmpFloat)

    resu = PID_Results( myeval.xyz, zeros(myeval.n), zeros(myeval.n), ones(myeval.n) )
    #                    support         q                ∇               p

    _lb = zeros(Float64,myeval.m)
    _ub = zeros(Float64,myeval.m)
    _l  = zeros(Float64,myeval.n)
    _u  = [Inf  for j=1:myeval.n]

    return PID_Data( 1/log(2), myeval, resu ,  pdummy, _lb, _ub, _l, _u)
    ;      #         bits                      _space
end

function pid!(pd::PID_Data, p::Vector{Float64} ) :: Void
    @assert length(p) == pd.myeval.n   "Length of p-vector doesn't match PID_Data"

    const n = pd.myeval.n
    const m = pd.myeval.m

    solver = Mosek.MosekSolver(MSK_IPAR_INTPNT_MULTI_THREAD=0, MSK_IPAR_INTPNT_MAX_ITERATIONS=500, MSK_DPAR_OPTIMIZER_MAX_TIME=100.0)

    const model = MathProgBase.NonlinearModel( solver )

    for i in 1:n
        x,y,z = pd.myeval.xyz[i]
        pd._space[ x,y,z ]  =  p[i]
    end

    set_rhs!(pd.myeval,pd._space)

    MathProgBase.loadproblem!(model, n, m, pd._l, pd._u, pd._lb, pd._ub, :Max, pd.myeval)

    MathProgBase.optimize!(model)

    if MathProgBase.status(model) ∈ [:Solve_Succeeded, :Optimal, :NearOptimal, :Suboptimal, :FeasibleApproximate]
        my_sol  = MathProgBase.getsolution(model)
        my_dual = MathProgBase.getconstrduals(model)

        for i in 1:n # changed m to n
            x,y,z = pd.results.support[i]
            pd.results.q[i]   =   my_sol[   pd.myeval.varidx[x,y,z]   ]
        end

        pd.results.λ .= 0.
        for k in 1:pd.myeval.m
            mr = pd.myeval.mr_eq[k]
            if mr[1] == "xy"
                x,y = mr[2], mr[3]
                for z in 1:pd.myeval.n_z
                    i = pd.myeval.varidx[x,y,z]
                    i==0  || (  pd.results.λ[i] += my_dual[k]  )
                end
            elseif mr[1]=="xz"
                x,z = mr[2], mr[3]
                for y in  1:pd.myeval.n_y
                    i = pd.myeval.varidx[x,y,z]
                    i==0  || (  pd.results.λ[i] += my_dual[k]  )
                end
            else
                @assert false
            end
        end
    else
        # Problematic: status == :UserLimit -- and of course all other kinds of crashes
        @show status(model)
        @assert status(model) ∈ [:Solve_Succeeded, :Optimal, :NearOptimal, :Suboptimal, :FeasibleApproximate]
    end
    ;
end #^ pid!()

end #module Opt_PID


text:
using ImageQuilting
using GeoStatsBase
using GeoStatsImages
using Statistics
using Plots, VisualRegressionTests
using Test, Pkg, Random

# workaround GR warnings
ENV["GKSwstype"] = "100"

# environment settings
islinux = Sys.islinux()
istravis = "TRAVIS" ∈ keys(ENV)
datadir = joinpath(@__DIR__,"data")
visualtests = !istravis || (istravis && islinux)
if !istravis
  Pkg.add("Gtk")
  using Gtk
end

@testset "ImageQuilting.jl" begin
  @testset "Basic checks" begin
    # the output of a homogeneous image is also homogeneous
    TI = ones(20,20,20)
    reals = iqsim(TI, (10,10,10), size(TI))
    @test reals[1] == TI

    # categories are obtained from training image only
    ncateg = 3; TI = rand(0:ncateg, 20, 20, 20)
    reals = iqsim(TI, (10,10,10), size(TI))
    @test Set(reals[1]) ⊆ Set(TI)
  end

  @testset "Soft data" begin
    # trends with soft data
    TI = [zeros(10,20,1); ones(10,20,1)]
    trend = [zeros(20,10,1) ones(20,10,1)]
    reals = iqsim(TI, (10,10,1), size(TI), soft=[(trend,TI)], tol=1)
    @test mean(reals[1][:,1:10,:]) ≤ mean(reals[1][:,11:20,:])

    # no side effects with soft data
    TI = ones(20,20,20)
    TI[:,5,:] .= NaN
    aux = fill(1.0, size(TI))
    iqsim(TI, (10,10,10), size(TI), soft=[(aux,aux)])
    @test aux == fill(1.0, size(TI))

    # auxiliary variable with integer type
    TI = ones(20,20,20)
    aux = [i for i in 1:20, j in 1:20, k in 1:20]
    iqsim(TI, (10,10,10), size(TI), soft=[(aux,aux)])
    @test aux == [i for i in 1:20, j in 1:20, k in 1:20]
  end

  @testset "Hard data" begin
    # hard data is honored everywhere
    TI = ones(20,20,20)
    obs = rand(size(TI)...)
    data = Dict(CartesianIndex(i,j,k)=>obs[i,j,k] for i=1:20, j=1:20, k=1:20)
    reals = iqsim(TI, (10,10,10), size(TI), hard=data)
    @test reals[1] == obs

    # multiple realizations with hard data
    TI = ones(20,20,20)
    data = Dict(CartesianIndex(20,20,20)=>10)
    reals = iqsim(TI, (10,10,10), size(TI), hard=data, nreal=3)
    for real in reals
      @test real[20,20,20] == 10
    end
  end

  @testset "Masked grids" begin
    # masked simulation domain
    TI = ones(20,20,20)
    shape = Dict{CartesianIndex{3},Real}()
    active = trues(size(TI))
    for i=1:20, j=1:20, k=1:20
      if (i-10)^2 + (j-10)^2 + (k-10)^2 < 25
        push!(shape, CartesianIndex(i,j,k)=>NaN)
        active[i,j,k] = false
      end
    end
    reals = iqsim(TI, (10,10,10), size(TI), hard=shape)
    @test all(isnan.(reals[1][.!active]))
    @test all(.!isnan.(reals[1][active]))

    # masked training image
    TI = ones(20,20,20)
    TI[:,5,:] .= NaN
    reals = iqsim(TI, (10,10,10), size(TI))
    @test reals[1] == fill(1.0, size(TI))
    TI[1,5,:] .= 0
    reals = iqsim(TI, (10,10,10), size(TI))
    @test reals[1] == fill(1.0, size(TI))

    # masked domain and masked training image
    TI = ones(20,20,20)
    TI[:,5,:] .= NaN
    aux = fill(1.0, size(TI))
    shape = Dict(CartesianIndex(i,j,k)=>NaN for i=1:20, j=5, k=1:20)
    reals = iqsim(TI, (10,10,10), size(TI), hard=shape)
    @test all(isnan.(reals[1][:,5,:]))
    @test all(reals[1][:,1:4,:] .== 1)
    @test all(reals[1][:,6:20,:] .== 1)
    reals = iqsim(TI, (10,10,10), size(TI), hard=shape, soft=[(aux,aux)])
    @test all(isnan.(reals[1][:,5,:]))
    @test all(reals[1][:,1:4,:] .== 1)
    @test all(reals[1][:,6:20,:] .== 1)
  end

  @testset "Minimum error cut" begin
    # 3D cut
    TI = ones(20,20,20)
    _, _, voxs = iqsim(TI, (10,10,10), overlap=(1/3,1/3,1/3), debug=true)
    @test 0 ≤ voxs[1] ≤ 1
  end

  @testset "Simulation paths" begin
    # different simulation paths
    for kind in [:raster,:dilation,:random]
      path = ImageQuilting.genpath((10,10,10), kind, Int[])
      @test length(path) == 1000
    end

    # data is visited first if present
    path = ImageQuilting.genpath((10,10,10), :data, [1,1000])
    @test path[1:2] == [1,1000] || path[1:2] == [1000,1]
  end

  @testset "Voxel reuse" begin
    # mean voxel reuse is in range [0,1]
    TI = rand(20,20,20)
    μ, σ = voxelreuse(TI, (10,10,10), nreal=1)
    @test 0 ≤ μ ≤ 1
  end

  if visualtests
    for TIname in ["Strebelle","StoneWall"]
      Random.seed!(2017)
      TI = training_image(TIname)[1:50,1:50,:]
      reals = iqsim(TI, (30,30,1), size(TI), nreal=4)
      ps = [heatmap(real[:,:,1]) for real in reals]
      @plottest plot(ps...) joinpath(datadir,"Reals"*TIname*".png") !istravis
    end
    for TIname in ["StoneWall","WalkerLake"]
      Random.seed!(2017)
      TI = training_image(TIname)[1:20,1:20,:]
      @plottest voxelreuseplot(TI) joinpath(datadir,"Voxel"*TIname*".png") !istravis
    end
  end

  @testset "GeoStats.jl API" begin
    geodata = PointSetData(Dict(:variable => [1.,0.,1.]), [25. 50. 75.; 25. 75. 50.])
    grid    = RegularGrid{Float64}(100,100)
    problem = SimulationProblem(geodata, grid, :variable, 3)

    TI = training_image("Strebelle")
    inactive = [CartesianIndex(i,j,1) for i in 1:30 for j in 1:30]
    solver = ImgQuilt(:variable => (TI=TI, tilesize=(30,30,1), inactive=inactive))

    Random.seed!(2017)
    solution = solve(problem, solver)
    @test keys(solution.realizations) ⊆ [:variable]

    incomplete_solver = ImgQuilt()
    @test_throws AssertionError solve(problem, incomplete_solver)

    if visualtests
      @plottest plot(solution,size=(1000,300)) joinpath(datadir,"GeoStatsAPI.png") !istravis
    end
  end

  # if ImageQuilting.cl ≠ nothing && ImageQuilting.clfft ≠ nothing
    # @testset "GPU support" begin
      # # CPU and GPU give same results
      # TI = ones(20,20,20)
      # TI[10:end,:,:] = 2
      # Random.seed!(0); realscpu = iqsim(TI, (10,10,10), size(TI), gpu=false)
      # Random.seed!(0); realsgpu = iqsim(TI, (10,10,10), size(TI), gpu=true)
      # @test realscpu[1] == realsgpu[1]
    # end
  # end
end


text:
# using PyPlot
include("nonidentifiable_logistic.jl")
include("asvar.jl")
include("zigzagtools.jl")

# data initialization
logistic_parameter = -1
dimension = 2
n_observations = 1000
stepsize_MALA = 0.012
dataX = randn(n_observations)
p_vec = 1./(1 + exp(-logistic_parameter * dataX))
dataY = convert(Vector{Int64},rand(n_observations) .<= p_vec)
rho = 1

# for asymptotic variance /ess
n_batches = 500

# MALA initialization
Xi_init = [-1.0;0.0]
n_iter_MALA = 1000
g(xi) = xi[1] + xi[2].^2
U_fun(xi) = rho/2*sumabs2(xi) + sum(log(1 + exp(g(xi).*dataX))- g(xi) * dataX .* dataY)
grad_U(xi) = [rho * xi[1] + sum((1./(1+exp(-g(xi).*dataX))-dataY) .* dataX);
              rho * xi[2] + 2 * xi[2] * sum((1./(1+exp(-g(xi).*dataX))-dataY).* dataX)]

# n_iter_MALA =  * n_iter_ZZ
tic()
Xi_MALA = MALA(U_fun, grad_U, stepsize_MALA, 2, n_iter_MALA, Xi_init)
t_mala = toc()
asvar_mala = (BatchMeans(Xi_MALA[1,:],n_batches=n_batches),BatchMeans(Xi_MALA[2,:],n_batches=n_batches))

# Zig-Zag initialization
tic()
(T,Xi_ZZ) = zzNonidentifiableLogistic(dataX, dataY, 100*n_iter_MALA, Xi_init, rho)
t_zz = toc();
# samples_zz = equidistantSamples(T,Xi_ZZ,n_iter_MALA)
# samples_zz_ext = equidistantSamples(T,Xi_ZZ, Int(1e7))
# asvar_zz = (BatchMeans(samples_zz[1,:],n_batches=n_batches), BatchMeans(samples_zz[2,:],n_batches=n_batches))
# ess_zz_extended = (ess(samples_zz_ext[1,:]),ess(samples_zz_ext[2,:]))

n_batches_consensus = 10
stepsize_consensus = stepsize_MALA * n_batches_consensus
Xi_consensus = consensusNonidentifiableLogistic(dataX, dataY, n_batches_consensus, stepsize_consensus, n_iter_MALA, Xi_init)

n_batches_SGLD = 10
stepsize_SGLD = 0.001
Xi_SGLD = SGLDNonidentifiableLogistic(dataX, dataY, n_batches_SGLD, stepsize_SGLD, n_iter_MALA)

# tic()
# (T_cv,Xi_cv) = @profile cvNonidentifiableLogistic(dataX, dataY, n_proposedswitches_CV, Xi_init, rho)
n_proposedswitches_CV = n_iter_MALA * 100
(T_cv,Xi_cv) = cvNonidentifiableLogistic(dataX, dataY, n_proposedswitches_CV, Xi_init, rho)
# t_cv = toc()
# samples_cv = equidistantSamples(T_cv,Xi_cv, n_iter_MALA)
# asvar_cv = (BatchMeans(samples_cv[1,:],n_batches=n_batches), BatchMeans(samples_cv[2,:],n_batches=n_batches))

# variance_estimate = (PooledVariance([asvar_mala[1], asvar_zz[1], asvar_cv[1]], [length(Xi_MALA[1,:]), length(samples_zz[1,:]), length(samples_cv[1,:])], [var(Xi_MALA[1,:]), var(samples_zz[1,:]), var(samples_cv[1,:])]),
  # PooledVariance([asvar_mala[2], asvar_zz[2], asvar_cv[2]], [length(Xi_MALA[2,:]), length(samples_zz[2,:]), length(samples_cv[2,:])], [var(Xi_MALA[2,:]), var(samples_zz[2,:]), var(samples_cv[2,:])]))
# ess_mala = (ess(Xi_MALA[1,:], variance = variance_estimate[1], n_batches=n_batches), ess(Xi_MALA[2,:], variance = variance_estimate[2], n_batches=n_batches))
# ess_zz = (ess(samples_zz[1,:],variance = variance_estimate[1], n_batches=n_batches),ess(samples_zz[2,:], variance = variance_estimate[2], n_batches=n_batches))
# ess_cv = (ess(samples_cv[1,:], variance = variance_estimate[1], n_batches=n_batches),ess(samples_cv[2,:], variance = variance_estimate[2], n_batches=n_batches))

# samples_cv_ext = equidistantSamples(T,Xi_cv, Int(1e7))
# ess_cv_extended = (ess(samples_cv_ext[1,:]),ess(samples_cv_ext[2,:]))
# println("--------------------")
# println("ESS MALA : ", ess_mala, ", runtime MALA : ", t_mala, "s")
# println("ESS ZZ   : ", ess_zz, ", runtime ZZ   : ", t_zz, "s")
# println("ESS CV   : ", ess_cv, ", runtime CV   : ", t_cv, "s")

# using Gadfly
# l1 = layer(x=Xi_ZZ[1,:], y=Xi_ZZ[2,:],Theme(default_color=color("red")),Geom.path,order=1)
# l2 = layer(x=Xi_MALA[1,:], y=Xi_MALA[2,:],Theme(default_color=color("blue")),Geom.point,order=2)
# p1 = plot(l1,l2,Coord.cartesian(aspect_ratio=1))
# draw(PDF("export.pdf", 10cm,10cm), p1)
n_gridpoints = 100
x = linspace(-3,0.5, n_gridpoints)
y = linspace(-1.5,1.5,n_gridpoints)

xgrid = repmat(x,1,n_gridpoints)
ygrid = repmat(y',n_gridpoints,1)
z = zeros(n_gridpoints, n_gridpoints)
for i in 1:n_gridpoints
  for j in 1:n_gridpoints
    z[i,j] = exp(-U_fun([x[i];y[j]]))
  end
end

using PyPlot
figure(1);
clf();
cp = contour(xgrid, ygrid, z, colors="black", linewidth=2.0)
PyPlot.plot(Xi_ZZ[1,:], Xi_ZZ[2,:],color="blue",zorder=6)
# PyPlot.scatter(samples_zz[1,:], samples_zz[2,:],color="blue",zorder=6)
axis([-3.0,0.5,-1.5,1.5])
axis(:equal)
# PyPlot.plot(Xi_cv[1,:], Xi_cv[2,:],color="black",zorder=3)
# PyPlot.scatter(samples_cv[1,:], samples_cv[2,:],color="black",zorder=3)

figure(2);clf();
cp = contour(xgrid, ygrid, z, colors="black", linewidth=2.0)
PyPlot.scatter(Xi_MALA[1,:], Xi_MALA[2,:], color="blue",zorder=5)
axis([-3.0,0.5,-1.5,1.5])
axis(:equal)

figure(3);clf();
cp = contour(xgrid, ygrid, z, colors="black", linewidth=2.0)
PyPlot.scatter(Xi_consensus[1,:], Xi_consensus[2,:], color="blue")
axis([-3.0,0.5,-1.5,1.5])
axis(:equal)

figure(4);clf();
cp = contour(xgrid, ygrid, z, colors="black", linewidth=2.0)
PyPlot.scatter(Xi_SGLD[1,:], Xi_SGLD[2,:], color="blue")
axis([-3.0,0.5,-1.5,1.5])
axis(:equal)

figure(5);clf();
cp = contour(xgrid, ygrid, z, colors="black", linewidth=2.0)
PyPlot.plot(Xi_cv[1,:], Xi_cv[2,:], color="blue")
axis([-3.0,0.5,-1.5,1.5])
axis(:equal)

figure(6);clf()
cp = contour(xgrid, ygrid, z, colors="black", linewidth=2.0)
axis([-3.0,0.5,-1.5,1.5])
axis(:equal)

max_stars_repo_name: jbierkens/zigzag-experiments

max_issues_repo_name: jbierkens/zigzag-experiments

max_forks_repo_name: jbierkens/zigzag-experiments


text:
using OEIS
using Base.Test

# prime numbers
s = oeis(:A000040)
@test s.values == primes(s.values[end])

# test Int128
s = oeis(:A000521)
@test typeof(s.values[16]) == Int128

# test handling of too many values
@test_throws ErrorException oeis("primes")

# test handling of negatives
s = oeis(:A000594)
@test s.values[2] == -24


text:
using SimpleDirectMediaLayer
const SDL2 = SimpleDirectMediaLayer
using Test

SDL2_pkg_dir = joinpath(@__DIR__, "..","..")
audio_example_assets_dir = joinpath(SDL2_pkg_dir, "src/examples/audio_example")

# check that an audio device if available
SDL2.Init(UInt32(SDL2.INIT_VIDEO))
device = SDL2.Mix_OpenAudio(Int32(22050), UInt16(SDL2.MIX_DEFAULT_FORMAT), Int32(2), Int32(1024) )

if device == 0
SDL2.Mix_CloseAudio()
SDL2.Quit()

@testset "Init+Quit" begin
# Test that you can init and quit SDL_mixer multiple times correctly.
@test 0 == SDL2.Init(UInt32(SDL2.INIT_VIDEO))
@test 0 == SDL2.Mix_OpenAudio(Int32(22050), UInt16(SDL2.MIX_DEFAULT_FORMAT), Int32(2), Int32(1024) )
SDL2.Mix_CloseAudio()
SDL2.Quit()

@test 0 == SDL2.Init(UInt32(SDL2.INIT_VIDEO))
@test 0 == SDL2.Mix_OpenAudio(Int32(22050), UInt16(SDL2.MIX_DEFAULT_FORMAT), Int32(2), Int32(1024) )
SDL2.Mix_CloseAudio()
SDL2.Quit()
end

@testset "Sounds" begin
SDL2.Init(UInt32(0))
SDL2.Mix_OpenAudio(Int32(22050), UInt16(SDL2.MIX_DEFAULT_FORMAT), Int32(2), Int32(1024) )

med = SDL2.Mix_LoadWAV( joinpath(audio_example_assets_dir,"medium.wav")  )
@test med != C_NULL
@test -1 != SDL2.Mix_PlayChannel( Int32(-1), med, Int32(0) )
# Test that can play multiple times successfully
# (Note that if a sound overlaps with a previous sound, it will play on a
# different channel. The return value is which channel it plays on.)
@test -1 != SDL2.Mix_PlayChannel( Int32(-1), med, Int32(0) )
@test -1 != SDL2.Mix_PlayChannel( Int32(-1), med, Int32(0) )

# Test different overlapping sounds
scratch = SDL2.Mix_LoadWAV( joinpath(audio_example_assets_dir,"scratch.wav")  )
@test scratch != C_NULL
@test -1 != SDL2.Mix_PlayChannel( Int32(-1), scratch, Int32(0) )
@test -1 != SDL2.Mix_PlayChannel( Int32(-1), med, Int32(0) )
@test -1 != SDL2.Mix_PlayChannel( Int32(-1), scratch, Int32(0) )

SDL2.Mix_CloseAudio()
SDL2.Quit()
end

@testset "Music" begin
SDL2.Init(UInt32(0))
SDL2.Mix_OpenAudio(Int32(22050), UInt16(SDL2.MIX_DEFAULT_FORMAT), Int32(2), Int32(1024) )

# Load the music
music = SDL2.Mix_LoadMUS( joinpath(audio_example_assets_dir,"beat.wav") );
@test music != C_NULL

# Test playing/pausing the music
@test 0 == SDL2.Mix_PlayMusic( music, Int32(-1) )
SDL2.Mix_PauseMusic()
SDL2.Mix_ResumeMusic()
@test 0 == SDL2.Mix_HaltMusic()

# Test noops if no music is playing.
@test 0 == SDL2.Mix_HaltMusic()
SDL2.Mix_PauseMusic()
SDL2.Mix_ResumeMusic()

# Test playing multiple times
@test 0 == SDL2.Mix_PlayMusic( music, Int32(-1) )
@test 0 == SDL2.Mix_PlayMusic( music, Int32(-1) )

SDL2.Mix_CloseAudio()
SDL2.Quit()
end

end
SDL2.Quit()

max_stars_repo_name: Kyjor/julia-game

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_name: Kyjor/julia-game

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_name: Kyjor/julia-game

max_forks_repo_licenses: MIT Unlicense


text:
"""
In this module we define a MockTensor which inherits from the NDTensors.TensorStorage.
This can be used as the store for ITensor objects, but only stores the dimensions and so
can be used for tracking the dimension of contractions which would not be possible if
storing all the data.
"""

using ITensors

export MockTensor

""" Tensor store struct that just tracks tensor dimensions"""
struct MockTensor{T,N} <: AbstractArray{T, N}
    size::NTuple{N, Int}
end

MockTensor(size::NTuple{N, Int64}) where N = MockTensor{ComplexF64, N}(size)

"""Overload functions from base to make MockTensor usable"""
Base.copy(a::MockTensor{T, N}) where {T, N} = MockTensor{T, N}(a.size)
Base.length(a::MockTensor) = prod(size(a))
Base.size(a::MockTensor) = a.size
Base.getindex(::MockTensor, ::Int64) = NaN
Base.getindex(::MockTensor, i...) = NaN
Base.show(io::IO, ::MIME"text/plain", a::MockTensor) = print(io, "MockTensor with dims $(Tuple(a.size))")