text:
####################################################
# HELLO KROME test evolves the network
#  FK1 -> FK2 -> FK3
# where FK* are fake species

# See the corresponding `test.f90` file for comparison of the implementation.
#
# `test_hello()` generates a file `julia.66`.
# The Julia reference values can be found in `julia.66.ref`.
# The Fortran reference values can be found in `fort.66.ref`.
#
# IMPORTANT: Many C bindings require to pass values as pointers. For example, `krome` cannot be
#            called with `Tgas` being of type `Float64`, but must be called with a **mutable**
#            value. Since scalars are immutable in Julia, the easiest is to create a
#            zero-dimensional array, which can be accomplished by calling `fill(<value>)`.

using KROME
using Printf

function test_hello(io::IO=stdout)
  krome_init() # init krome (mandatory)

  nmols = krome_nmols()[] # read Fortran module variable
  x = zeros(nmols) # default abundances (number density)
  idx_FK1 = krome_idx_FK1()[] + 1 # Julia has 1-based indices
  x[idx_FK1] = 1.0 # FK1 initial abundance (number density)

  Tgas = fill(1e2) # gas temperature, not used (K)
  dt = fill(1e-5) # time-step (arbitrary)
  t = 0.0 # time (arbitrary)

  println(io, "#time  FK1 FK2 FK3")
  while true
    dt[] = dt[] * 1.1 # increase timestep
    krome(x, Tgas, dt) # call KROME
    print(io, fsn(t))
    for value in x
      print(io, fsn(value))
    end
    println(io)
    t = t + dt[] # increase time
    t > 5 && break # break loop after 5 time units
  end
end

# Return formatted value in Fortran scientific notation
function fsn(value, width=17, precision=8)
  @assert precision <= 64 "precision must be <= 64"

  # Extract exponent and base
  if value == 0.0
    exponent = -1
    base = 0.0
  else
    exponent = floor(Int, log10(abs(value)))
    base = value / float(10)^exponent
  end

  # Make base have a leading zero
  base /= 10
  exponent += 1

  # Get formatted values
  base_formatted = (@sprintf "%.64f" base)[1:(precision + 2 + (base < 0))]
  exponent_formatted = @sprintf "E%+04d" exponent

  padding = width - length(base_formatted) - length(exponent_formatted)

  return " "^padding * base_formatted * exponent_formatted
end


text:
struct VTPRUNIFACParam <: EoSParam
    A::PairParam{Float64}
    B::PairParam{Float64}
    C::PairParam{Float64}
    Q::SingleParam{Float64}
end

abstract type VTPRUNIFACModel <: UNIFACModel end

struct VTPRUNIFAC{c<:EoSModel} <: VTPRUNIFACModel
    components::Array{String,1}
    icomponents::UnitRange{Int}
    groups::GroupParam
    params::VTPRUNIFACParam
    puremodel::Vector{c}
    absolutetolerance::Float64
    references::Array{String,1}
end

@registermodel VTPRUNIFAC
export VTPRUNIFAC

function VTPRUNIFAC(components; puremodel=PR,
    userlocations=String[],
     verbose=false)
    groups = GroupParam(components, ["Activity/UNIFAC/VTPR/VTPR_groups.csv"]; verbose=verbose)

    params = getparams(groups, ["Activity/UNIFAC/VTPR/VTPR_like.csv", "Activity/UNIFAC/VTPR/VTPR_unlike.csv"]; userlocations=userlocations,  asymmetricparams=["A","B","C"], ignore_missing_singleparams=["A","B","C"], verbose=verbose)
    A  = params["A"]
    B  = params["B"]
    C  = params["C"]
    Q  = params["Q"]
    icomponents = 1:length(components)

    init_puremodel = [puremodel([groups.components[i]]) for i in icomponents]
    packagedparams = VTPRUNIFACParam(A,B,C,Q)
    references = String[]
    model = VTPRUNIFAC(components,icomponents,groups,packagedparams,init_puremodel,1e-12,references)
    return model
end

function activity_coefficient(model::VTPRUNIFACModel,V,T,z)
    return exp.(@f(lnγ_res))
end

ext: jl


text:
module Merge
import Base: start, next, done, eltype, length
using DataStructures
export merge_sorted

immutable SortedMerger{I}
    keyfunc::Function
    its::Vector{I}
end

function merge_sorted(its, f=identity)
    SortedMerger(f, its)
end

function start(it::SortedMerger)
    el_type = eltype(eltype(it.its))
    if el_type <: AbstractString
        el_type = AbstractString
    end
    h = binary_minheap(Tuple{el_type, Int})
    l = length(it.its)
    its = Vector(l)
    values = Vector{eltype(eltype(it.its))}(l)
    for (i, xs) in enumerate(it.its)
        s = start(xs)
        if done(xs, s)
            continue
        end

        (x, s) = next(xs, s)
        key = it.keyfunc(x)
        values[i] = x
        its[i] = (xs, s)
        push!(h, (key, i))
    end
    (h, its, values)
end

function next(it::SortedMerger, state)
    (h, its, values) = state
    (key, i) = pop!(h)
    next_x = values[i]
    (xs, s) = its[i]
    if !done(xs, s)
        (x, s) = next(xs, s)
        key = it.keyfunc(x)
        its[i] = (xs, s)
        values[i] = x
        push!(h, (key, i))
    end
    (next_x, (h, its, values))
end

function done(it::SortedMerger, state)
    (h, its, values) = state
    isempty(h)
end

end # module


text:
"""
    ImpactContact{T,N} <: Contact{T,N}

    contact object for impact (i.e., no friction)

    collision: Collision
"""
mutable struct ImpactContact{T,N} <: Contact{T,N}
    friction_parameterization::SMatrix{0,2,T,0}
    collision::Collision{T,0,3,0}
end

# function Base.show(io::IO, mime::MIME{Symbol("text/plain")}, contact::ImpactContact)
#     summary(io, contact)
#     println(io, "")
#     println(io, "friction_parameterization: "*string(contact.friction_parameterization))
#     println(io, "collision:                 "*string(contact.collision))
# end

function ImpactContact(body::Body{T}, normal::AbstractVector{T};
    contact_origin=szeros(T, 3),
    contact_radius=0.0) where T

    # contact directions
    V1, V2, V3 = orthogonal_columns(normal) #
    A = [V1 V2 V3]
    Ainv = inv(A)
    contact_normal = Ainv[3, SA[1; 2; 3]]'

    # friction parametrization
    parameterization = szeros(T, 0, 2)

    # collision
    collision = SphereHalfSpaceCollision(szeros(T, 0, 3), contact_normal, SVector{3}(contact_origin), contact_radius)

    ImpactContact{Float64,2}(parameterization, collision)
end

function constraint(mechanism, contact::ContactConstraint{T,N,Nc,Cs}) where {T,N,Nc,Cs<:ImpactContact{T,N}}
    # contact model
    model = contact.model

    # parent
    pbody = get_body(mechanism, contact.parent_id)
    xp, qp = next_configuration(pbody.state, mechanism.timestep)

    # child
    cbody = get_body(mechanism, contact.child_id)
    xc, qc = next_configuration(cbody.state, mechanism.timestep)

    # constraint
    SVector{1,T}(distance(model.collision, xp, qp, xc, qc) - contact.impulses_dual[2][1])
end

function constraint_jacobian(contact::ContactConstraint{T,N,Nc,Cs,N½}) where {T,N,Nc,Cs<:ImpactContact{T,N},N½}
    γ = contact.impulses[2] + REG * neutral_vector(contact.model)
    s = contact.impulses_dual[2] + REG * neutral_vector(contact.model)
    ∇s = hcat(γ, -Diagonal(sones(N½)))
    ∇γ = hcat(s, -Diagonal(szeros(N½)))
    return [∇s ∇γ]
end

function constraint_jacobian_configuration(relative::Symbol, model::ImpactContact,
    xp::AbstractVector, vp::AbstractVector, qp::Quaternion, ωp::AbstractVector,
    xc::AbstractVector, vc::AbstractVector, qc::Quaternion, ωc::AbstractVector,
    timestep)

    X = ∂distance∂x(relative, model.collision, xp, qp, xc, qc)
    Q = ∂distance∂q(relative, model.collision, xp, qp, xc, qc)

    return [X Q]
end

function constraint_jacobian_velocity(relative::Symbol, model::ImpactContact,
    xp::AbstractVector, vp::AbstractVector, qp::Quaternion, ωp::AbstractVector,
    xc::AbstractVector, vc::AbstractVector, qc::Quaternion, ωc::AbstractVector,
    timestep)

    # recover current orientation
    if relative == :parent
        x = next_position(xp, -vp, timestep)
        ∂x∂v = linear_integrator_jacobian_velocity(x, vp, timestep)
        q = next_orientation(qp, -ωp, timestep)
        ∂q∂ω = rotational_integrator_jacobian_velocity(q, ωp, timestep)
    elseif relative == :child
        x = next_position(xc, -vc, timestep)
        ∂x∂v = linear_integrator_jacobian_velocity(x, vc, timestep)
        q = next_orientation(qc, -ωc, timestep)
        ∂q∂ω = rotational_integrator_jacobian_velocity(q, ωc, timestep)
    end

    # Jacobian
    V = ∂distance∂x(relative, model.collision, xp, qp, xc, qc) * ∂x∂v
    Ω = ∂distance∂q(relative, model.collision, xp, qp, xc, qc) * ∂q∂ω

    return [V Ω]
end

function force_mapping(relative::Symbol, model::ImpactContact,
    xp::AbstractVector, qp::Quaternion,
    xc::AbstractVector, qc::Quaternion)

    X = contact_normal(model.collision, xp, qp, xc, qc)'

    if relative == :parent
        return X
    elseif relative == :child
        return -1.0 * X
    end
end

function ∂force_mapping_jvp∂x(relative::Symbol, jacobian::Symbol,
    model::ImpactContact,
    xp::AbstractVector, qp::Quaternion,
    xc::AbstractVector, qc::Quaternion,
    λ::AbstractVector)

    X = ∂contact_normal_transpose∂x(jacobian, model.collision, xp, qp, xc, qc) * λ[1]

    if relative == :parent
        return X
    elseif relative == :child
        return -1.0 * X
    end
end

function ∂force_mapping_jvp∂q(relative::Symbol, jacobian::Symbol,
    model::ImpactContact,
    xp::AbstractVector, qp::Quaternion,
    xc::AbstractVector, qc::Quaternion,
    λ::AbstractVector)

    X = ∂contact_normal_transpose∂q(jacobian, model.collision, xp, qp, xc, qc) * λ[1]

    if relative == :parent
        return X
    elseif relative == :child
        return -1.0 * X
    end
end


text:
using GraphQLClient

@testset "query" begin
    gateway_url = "https://gateway.thegraph.com/network"

    @testset "verify_ipfshashes" begin
        # Should fail due to bad prefix
        hashes = [
            "QmauYgPmss6CEZXtaRvvGW2oiyLqxpoCkWNCTmFPVTFDfk",
            "AmauYgPmss6CEZXtaRvvGW2oiyLqxpoCkWNCTmFPVTFDfk",
        ]
        @test !verify_ipfshashes(hashes)

        # Should pass for good hashes
        hashes = [
            "QmauYgPmss6CEZXtaRvvGW2oiyLqxpoCkWNCTmFPVTFDfk",
            "QmhiYgPmss6CEZXtaRvvGW2oiyLqxpoCkWNCTmFPVTFDfk",
        ]
        @test verify_ipfshashes(hashes)

        # Should fail for empty strings
        hashes = [
            "QmauYgPmss6CEZXtaRvvGW2oiyLqxpoCkWNCTmFPVTFDfk",
            "QmhiYgPmss6CEZXtaRvvGW2oiyLqxpoCkWNCTmFPVTFDfk",
            "",
        ]
        @test !verify_ipfshashes(hashes)

        # Should pass for empty lists
        hashes = String[]
        @test verify_ipfshashes(hashes)
    end

    @testset "ipfshash_in" begin
        # Should handle nominal input
        whitelist = ["a", "b", "c"]
        pinnedlist = ["d"]
        output = ipfshash_in(whitelist, pinnedlist)
        @test output == ["a", "b", "c", "d"]

        # Should handle empty inputs
        whitelist = String[]
        pinnedlist = String[]
        output = ipfshash_in(whitelist, pinnedlist)
        @test output == []

        # Should handle overlapping input
        whitelist = ["a", "b", "c"]
        pinnedlist = ["c", "d"]
        output = ipfshash_in(whitelist, pinnedlist)
        @test output == ["a", "b", "c", "d"]
    end

    @testset "ipfshash_not_in" begin
        # Should handle nominal input
        blacklist = ["a", "b", "c"]
        frozenlist = ["d"]
        output = ipfshash_not_in(blacklist, frozenlist)
        @test output == ["a", "b", "c", "d"]

        # Should handle empty inputs
        blacklist = String[]
        frozenlist = String[]
        output = ipfshash_not_in(blacklist, frozenlist)
        @test output == []

        # Should handle overlapping input
        blacklist = ["a", "b", "c"]
        frozenlist = ["c", "d"]
        output = ipfshash_not_in(blacklist, frozenlist)
        @test output == ["a", "b", "c", "d"]
    end

    @testset "snapshot ipfs filtering subgraph deployments" begin
        client = Client(gateway_url)
        all_hashes = [
            deployment["ipfsHash"] for deployment in
            query(client, "subgraphDeployments"; query_args=Dict("first" => 1000, "where" => Dict("signalledTokens_gte" => "1000000000000000000000")), output_fields="ipfsHash").data["subgraphDeployments"]
        ]

        # Should return the one whitelisted subgraph
        repo = snapshot(client, [all_hashes[1]], String[])
        @test length(repo.subgraphs) == 1
        @test repo.subgraphs[1].ipfshash == all_hashes[1]

        # Should return the one subgraph not blacklisted
        repo = snapshot(client, String[], all_hashes[2:end])
        @test length(repo.subgraphs) == 1
        @test repo.subgraphs[1].ipfshash == all_hashes[1]

        # Should return all subgraphs
        repo = snapshot(client, String[], String[])
        @test length(repo.subgraphs) == length(all_hashes)

        # Should return subgraphs in whitelist that aren't in blacklist
        repo = snapshot(client, all_hashes, all_hashes[2:end])
        @test length(repo.subgraphs) == 1
        @test repo.subgraphs[1].ipfshash == all_hashes[1]
    end

    @testset "snapshot ipfs filtering indexer allocations" begin
        client = Client(gateway_url)
        all_hashes = [
            deployment["ipfsHash"] for deployment in
            query(client, "subgraphDeployments"; query_args=Dict("first" => 1000, "where" => Dict("signalledTokens_gte" => "1000000000000000000000")), output_fields="ipfsHash").data["subgraphDeployments"]
        ]

        # Should have only one possible subgraph to allocate to
        repo = snapshot(client, [all_hashes[1]], String[])
        all_allocations = unique([
            a -> a.ipfshash for indexer in repo.indexers for a in indexer.allocations
        ])
        @test length(all_allocations) == 1

        # Should have no allocations
        repo = snapshot(client, String[], all_hashes)
        @test sum(map(x -> length(x.allocations), repo.indexers)) == 0
    end

    @testset "frozen stake" begin
        # calculate the stake constraint for indexer
        client = Client(gateway_url)
        subgraphs = query_subgraphs(client, String[], String[])
        indexers = query_indexers(client, subgraphs)
        indexer = indexers[findfirst(i -> length(i.allocations) > 0, indexers)]
        id = indexer.id
        freezed_allocation = indexer.allocations[1]

        @test frozen_stake(client, id, String[]) == 0.0
        @test (frozen_stake(client, id, String[])) <
            frozen_stake(client, id, [freezed_allocation.ipfshash])
    end

    @testset "query indexer allocations" begin
        client = Client(gateway_url)
        subgraphs = query_subgraphs(client, String[], String[])
        indexers = query_indexers(client, subgraphs)
        indexer = indexers[findfirst(i -> length(i.allocations) > 0, indexers)]
        alloc_length = length(indexer.allocations)
        id = indexer.id

        # Should return the correct number of allocations
        allocations = query_indexer_allocations(client, id)
        @test length(allocations) == alloc_length
    end
end


text:
using MshReader,
      WriteVTK


text:
#http://pages.cs.wisc.edu/~psilord/blog/data/chess-pages/nonsliding.html
const K_CLEAR_FILES = (0x7f7f7f7f7f7f7f7f, 0xffffffffffffffff,
                       0xfefefefefefefefe, 0xfefefefefefefefe,
                       0xfefefefefefefefe, 0xffffffffffffffff,
                       0x7f7f7f7f7f7f7f7f, 0x7f7f7f7f7f7f7f7f)
const K_MASK_FILES = (0x8080808080808080, 0xffffffffffffffff,
                      0x0101010101010101, 0x0101010101010101,
                      0x0101010101010101, 0xffffffffffffffff,
                      0x8080808080808080, 0x8080808080808080)
const K_SHIFTS = (9, 8, 7, -1, -9, -8, -7, 1)

#pseudo moves
function kingMovesGen()
    king_moves = Dict{UInt64,UInt64}()

    for square in values(PGN2UINT)
        targets = EMPTY
        for (mask, clear, shift) in zip(K_MASK_FILES, K_CLEAR_FILES, K_SHIFTS)
            if square & mask != EMPTY
                targets |= (square & clear) << shift
            else
                targets |= square << shift
            end
        end
        targets != EMPTY ? push!(king_moves, square=>targets) : continue
    end

    return king_moves
end
const KING = kingMovesGen()

function inCheck(b::Board, white::Bool, K::UInt64=EMPTY)
    if K == EMPTY; white ? K = b.white.K : K = b.black.K end
    white ? enemy = b.black : enemy = b.white

    if KNIGHT[K] & enemy.N != EMPTY; return true end

    if KING[K] & enemy.K != EMPTY; return true end

    ortho = getMagicAttack(K, b.taken, true)
    if ortho & enemy.R != EMPTY || ortho & enemy.Q != EMPTY; return true end

    diago = getMagicAttack(K, b.taken, false)
    if diago & enemy.B != EMPTY || diago & enemy.Q != EMPTY; return true end

    white ? x_pawn = PAWN_X_WHITE : x_pawn = PAWN_X_BLACK
    if x_pawn[K] & enemy.P != EMPTY; return true end

    return false
end

function getCastlingMoves!(moves::Moves, K::UInt64, b::Board, white::Bool)
    if b.castling == NOCASTLING || inCheck(b, white); return end
    if white
        if CQ&b.castling != NOCASTLING && A1&b.white.R != EMPTY
            if B1C1D1 & b.taken == EMPTY
                if ~inCheck(b, white, D1) && ~inCheck(b, white, C1)
                    push!(moves, Move(:king, K, C1, NONE, EMPTY, :none, CQ))
                end
            end
        end
        if CK&b.castling != NOCASTLING && H1&b.white.R != EMPTY
            if F1G1 & b.taken == EMPTY
                if ~inCheck(b, white, F1) && ~inCheck(b, white, G1)
                    push!(moves, Move(:king, K, G1, NONE, EMPTY, :none, CK))
                end
            end
        end
    else
        if Cq&b.castling != NOCASTLING && A8&b.black.R != EMPTY
            if B8C8D8 & b.taken == EMPTY
                if ~inCheck(b, white, D8) && ~inCheck(b, white, C8)
                    push!(moves, Move(:king, K, C8, NONE, EMPTY, :none, Cq))
                end
            end
        end
        if Ck&b.castling != NOCASTLING && H8&b.black.R != EMPTY
            if F8G8 & b.taken == EMPTY
                if ~inCheck(b, white, F8) && ~inCheck(b, white, G8)
                    push!(moves, Move(:king, K, G8, NONE, EMPTY, :none, Ck))
                end
            end
        end
    end
end


text:
module CollabUtils2SilverBlock

using ..Ahorn, Maple

@mapdef Entity "CollabUtils2/SilverBlock" SilverBlock(x::Integer, y::Integer, width::Integer=Maple.defaultBlockWidth, height::Integer=Maple.defaultBlockHeight)

const placements = Ahorn.PlacementDict(
    "Silver Block (Collab Utils 2)" => Ahorn.EntityPlacement(
        SilverBlock,
        "rectangle"
    ),
)

frame = "CollabUtils2/silverblock"
strawberry = "CollabUtils2/silverBerry/idle00"

Ahorn.minimumSize(entity::SilverBlock) = 16, 16
Ahorn.resizable(entity::SilverBlock) = true, true

Ahorn.selection(entity::SilverBlock) = Ahorn.getEntityRectangle(entity)

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::SilverBlock, room::Maple.Room)
    x, y = Ahorn.position(entity)

    width = Int(get(entity.data, "width", 32))
    height = Int(get(entity.data, "height", 32))

    tilesWidth = div(width, 8)
    tilesHeight = div(height, 8)

    for i in 2:tilesWidth - 1
        Ahorn.drawImage(ctx, frame, (i - 1) * 8, 0, 8, 0, 8, 8)
        Ahorn.drawImage(ctx, frame, (i - 1) * 8, height - 8, 8, 16, 8, 8)
    end

    for i in 2:tilesHeight - 1
        Ahorn.drawImage(ctx, frame, 0, (i - 1) * 8, 0, 8, 8, 8)
        Ahorn.drawImage(ctx, frame, width - 8, (i - 1) * 8, 16, 8, 8, 8)
    end

    for i in 2:tilesWidth - 1, j in 2:tilesHeight - 1
        Ahorn.drawImage(ctx, frame, (i - 1) * 8, (j - 1) * 8, 8, 8, 8, 8)
    end

    Ahorn.drawImage(ctx, frame, 0, 0, 0, 0, 8, 8)
    Ahorn.drawImage(ctx, frame, width - 8, 0, 16, 0, 8, 8)
    Ahorn.drawImage(ctx, frame, 0, height - 8, 0, 16, 8, 8)
    Ahorn.drawImage(ctx, frame, width - 8, height - 8, 16, 16, 8, 8)

    Ahorn.drawSprite(ctx, strawberry, width / 2, height / 2)
end

end


text:
using jInv.Mesh;
using MAT;
using jInv
using jInv.InverseSolve
using ShapeReconstructionPaLS
using ParamLevelSet;
using ShapeReconstructionPaLS.Utils;
using ShapeReconstructionPaLS.PointCloud;
using DelimitedFiles
using LinearAlgebra
using Statistics
using Distributed
using SparseArrays
using Random

ENV["MPLBACKEND"] = "Qt4Agg"

function ind2subv(shape, indices)
    CI = CartesianIndices(shape)
    return getindex.(Ref(CI), indices)
end

plotting = false;

if plotting
	using PyPlot;
	using jInvVisPyPlot
	close("all");
end

invertDip = true;
invertVis = false;

locateRBFwithGrads = true;
n = [128,128,128];
mesh = getRegularMesh([0.0 1.5 0.0 1.5 0.0 1.5],n);
n_tup = collect(n);
n_tup = tuple(n_tup...)

################################################
### Reading the model and pad ##################
################################################

model = "fandisksmall";
file = open(string(pwd(),"/../models/",model,".xyz"));

#file = open("pc_full.xyz")
A = readlines(file);
A = split.(A);
B =  Array{Array{Float64,1},1}(undef,length(A));
N = Array{Array{Float64,1},1}(undef,length(A));
for i = 1:1:size(A,1)
	 global B[i] = parse.(Float64,A[i])[1:3] .+ 0.75;
	 global N[i] = parse.(Float64,A[i])[4:6];
end
close(file);
A = Array{Float64,2}(transpose(hcat(B...)));
Normals = Array{Float64,2}(transpose(hcat(N...)));

##Sort the PC:
PC  = [ A Normals];
PC = PC[sortperm(PC[:, 1]), :];
writedlm("SortedPC.txt",PC);
println("PC:",size(PC)," ", PC[end,:])

npc = size(A,1); #Number of points in point cloud
B = [];

P = A;
println("Size of P:",size(P));

writedlm(string("originalPC.dat"),convert(Array{Float16},A));

#Clear vars:
indices = [];
A = []; tmp = [];
linIdx = [];

noiseAnglesDeg = 1;
#original was mean(2*Mesh.h);
noiseTrans = 0.015;
noiseAngles = deg2rad(noiseAnglesDeg);
#########################################################################################################
### Dipping Data ########################################################################################
#########################################################################################################

method = RBF10Based; methodName="RBF10Based";

nWorkers = nworkers();
println("nworkers=",nWorkers)
misfun = PCFun; ## least squares

	dipDataFilename = string("pointCloudData",model,"_dataRes",n,1);
	trans = [0.0 0 0 ; 0.0 0 0 ]
	theta_phi_dip = [0.0 0.0 ; 0.0 0.0];
	noiseTrans = 0.0;
	noiseAngles = deg2rad(0.0);
	npc = 2;
	prepareSyntheticPointCloudData(PC,mesh,npc,theta_phi_dip,trans,noiseTrans,noiseAngles,dipDataFilename);

	#########################################################################################################
	### Read the dipping data ###############################################################################
	#########################################################################################################
	Data,P,Normals,trans,theta_phi_dip,npc = readPointCloudDataFile(dipDataFilename);
	### Set up the dip inversion
	P =  Array{Array{Float64,2}}(P);
	Normals = Array{Array{Float64,2}}(Normals);
	pForDip = getPointCloudParam(mesh,P,Normals,theta_phi_dip,trans,nWorkers,method);
	ndips = npc;
	if(locateRBFwithGrads)
		pForDip_MATFree = getPointCloudParam(mesh,P,Normals,theta_phi_dip,trans,nWorkers,MATFree);
	end

	### Create Dip param if we use the GRAD of MATFREE to locate new RBFs
	dobsDirect = divideDataToWorkersData(nWorkers,Data);
	println("size of dobs",size(dobsDirect))
	#nnz = filter(x -> x != 0,dobsDirect[1]);
	#println("length of dobsdirect: ",length(dobsDirect[1]),"_",length(nnz));
	Wd_Dip = Array{Array{Float32}}(undef,nWorkers);
	for k=1:nWorkers
		Wd_Dip[k] = ones(Float32,size(dobsDirect[k]));
	end
	println("Wd dip:",size(Wd_Dip))

n_Moves_all = npc;
println("npc = ",npc);

isSimple = 0;
isRBF10 = 0;

###############################################################################################
############## USE RBF dictionary #############################################################
###############################################################################################

	nRBF = 2;
	isRBF10 = (method == RBF10BasedSimple1 || method==RBF10BasedSimple2 || method==RBF10Based);
	isSimple = (method==RBFBasedSimple1 || method==RBF10BasedSimple1 || method==RBFBasedSimple2 || method==RBF10BasedSimple2);

	if isRBF10
		numParamOfRBF = 10;
	else
		numParamOfRBF = 5;
	end
	n_m_simple = numParamOfRBF*nRBF;
	nAll = isSimple ? numParamOfRBF*nRBF : numParamOfRBF*nRBF + n_Moves_all*5;

	modfun = identityMod;

	if isSimple
		modfunForPlotting = (m)->MeshFreeParamLevelSetModelFunc(mesh,m ,computeJacobian = 0,
			sigma = getDefaultHeaviside() ,bf = 1, numParamOfRBF = numParamOfRBF);
	else
		modfunForPlotting = (m)->MeshFreeParamLevelSetModelFunc(mesh,m[1:(numParamOfRBF*nRBF)],computeJacobian = 0,
			sigma = getDefaultHeaviside() ,bf = 1, numParamOfRBF = numParamOfRBF);
	end

	Iact = 1.0;
	IactPlot = 1.0;
	sback = zero(Float32);
	a = -100000000.0;
	b = 10000000.0;
	boundsHigh = zeros(Float32,nAll) .+ b;
	# boundsHigh[end-5*npc+1:end] .= 0.01;

	boundsLow = zeros(Float32,nAll) .+ a;
	# boundsLow[end-5*npc+1:end] .= -0.01;
	# m0 - initial guess. mref - reference for regularization.
	mref = zeros(Float32,nAll);
	m0   = zeros(Float64,n_m_simple);
	mref[1:numParamOfRBF:n_m_simple] .= ParamLevelSet.centerHeavySide*2.0;
	if isRBF10
		# here we initialize B as an identity matrix
		mref[2:10:n_m_simple] .= 15.0;
		mref[5:10:n_m_simple] .= 15.0;
		mref[7:10:n_m_simple] .= 15.0;
		boundsLow[2:10:n_m_simple] .= 0.05;
		boundsLow[5:10:n_m_simple] .= 0.05;
		boundsLow[7:10:n_m_simple] .= 0.05;
		#boundsLow[end-9:end] .= -0.01
		#boundsHigh[end-9:end] .= 0.01;
	else
		mref[2:5:n_m_simple] .= 4.0;
		boundsLow[2:5:n_m_simple] .= 0.05;
	end
	mref[(numParamOfRBF-2):numParamOfRBF:n_m_simple] .= 0.75;
	mref[(numParamOfRBF-1):numParamOfRBF:n_m_simple] .= 0.75;
	mref[numParamOfRBF:numParamOfRBF:n_m_simple] .= 0.75;
	m0[:] .= mref[1:n_m_simple];

	if !isSimple
		m0 = wrapRBFparamAndRotationsTranslations(m0,theta_phi_dip,trans);
	end
	mref[:] = m0;
	mref[end-5*npc+1:end].=0.0;
	m0[(numParamOfRBF-2):numParamOfRBF:n_m_simple] .+= 0.1*randn(nRBF);
	m0[(numParamOfRBF-1):numParamOfRBF:n_m_simple] .+= 0.1*randn(nRBF);
	m0[numParamOfRBF:numParamOfRBF:n_m_simple]     .+= 0.1*randn(nRBF);

	alpha = 5e-4;

	II = (sparse(1.0I, nAll,nAll));
	regfun = (m, mref, M)->TikhonovReg(m,mref,M,II);
	HesPrec = getEmptyRegularizationPreconditioner();

	if isRBF10
		spd_reg = (m,mref,M) -> RBF_SPD_regularization(m,mref,nRBF);
		regfun 	 		= [regfun;spd_reg];
		alpha 		 	= [alpha;1e-1*alpha]; #was 1e-4 for cube\prism !!!!
		mref 			= [mref zeros(size(mref))];
	end

pMisRFs = getMisfitParam(pForDip, Wd_Dip, dobsDirect, misfun, Iact,sback);

cgit  = 100; # not used.
maxit = 10;
pcgTol = 1e-2; #not used

function myDump(mc::Vector,Dc,iter,pInv,pMis,method="",ndips=0,noiseAngle=0,noiseTranslation=0,useVisual=false,plotFun=modfunForPlotting)
	# if resultsFilename!=""
		# Temp = splitext(resultsFilename);
		# Temp = string(Temp[1],"_GN",iter,Temp[2]);

			#ntup = tuple((pInv.MInv.n)...);
		ntup = tuple(n...)
		fullMc = plotFun(mc)[1];
		#println("size of fullmc:",size(fullMc))
		writedlm(string("FandiskDirectParams_iter_",iter,".dat"),convert(Array{Float64},mc));
		writedlm(string("FandiskDirectVolume_iter_",iter,"_",ntup,".dat"),convert(Array{Float16},fullMc));
		if plotting
			fullMc = reshape(convert(Array{Float32},fullMc[:]),ntup);
			close(888);
			figure(888);
			plotModel(fullMc,includeMeshInfo = true,M_regular = pInv.MInv);
			savefig(string("FandiskDirectIt",iter,"_",ntup,".png"));
			sleep(0.5);
			# writedlm(string("FandiskDirectParams_iter_",iter,".dat"),convert(Array{Float64},mc));
			# writedlm(string("FandiskDirectVolume_iter_",iter,"_",ntup,".dat"),convert(Array{Float16},fullMc));

		end
	# end
end
pInv = getInverseParam(mesh,modfun,regfun,alpha,mref,copy(boundsLow),copy(boundsHigh),
                     maxStep = 2e-1,pcgMaxIter=cgit,pcgTol=pcgTol,
					 minUpdate = 0.01, maxIter = maxit,HesPrec=HesPrec);

# end

#### Projected Gauss Newton
mc = m0;
pInv.maxIter = 100;
Dc = nothing;
mc,Dc,flag,his = projGN(mc,pInv,pMisRFs,solveGN=projGNexplicit);
#mc[end-5*npc+1:end] .= 0.0;
print("DC shape:",size(Dc))

pInv.maxIter = 50;
myDump(mc,0,1,pInv,0,"",0,noiseAnglesDeg,noiseTrans,invertVis);

nRBF=nRBF;
times=0;

if !(method == MATBased || method == MATFree)
	outerIter = 80;
	grad_u = nothing;
	new_RBF_location = [];
	@sync for iterNum=2:outerIter
		println("Inversion using ", nRBF," basis functions")
		new_nRBF = 5;
		if (times > 100)
		  #times=0;
		  new_nRBF = 0;
		end
		global times = times + 1;

		# Dc - the data using our method.
		global Dc = Dc;

		#Lets compute the gradient of the misfit function:
		locateRBFwithGrads = true;
		print(mc[end-5*npc+1:end])
		if(locateRBFwithGrads)
			#(m11,theta_phi1,b1) 	= splitRBFparamAndRotationsTranslations(mc,nRBF,ndips,numParamOfRBF)
			#updateAngles(pForDip_MATFree,theta_phi1);
			pMis_Free = getMisfitParam(pForDip_MATFree, Wd_Dip, dobsDirect, misfun, Iact,sback);
			sigmaH = getDefaultHeaviside();
			allPoints = [P[1] ; P[2]];
			Dd, = computeMisfit(mc,pMisRFs);
			gradMis = (fetch(Dd[1])-Data).^2;
			sp = sortperm(gradMis,rev = true);
			sp = sp[1:round(Int64,length(sp)/4.0)];
			sp = sp[randperm(length(sp))];
			idxs = sp[1:new_nRBF];
			bwdSubs = [P[1] ; P[2]];
			println("all P size:",size(allPoints));
			Ii = allPoints[idxs,:];
		else
			u = modfunForPlotting(mc)[1];
			#Ii = findall(x -> x >0.2 && x<0.5,u);
			Ii = findall(x -> x >0.4 && x<0.8,u);
		end

		if(locateRBFwithGrads == false)
		if length(Ii) > new_nRBF
		Ii = Ii[randperm(length(Ii))[1:new_nRBF]];
		end
		end

    #while(Ii in new_RBF_location || Ii.-1 in new_RBF_location || Ii.+1 in new_RBF_location)
    #  Ii = findall(x -> x .>0.2 && x.<0.5,u);
    #  Ii = Ii[randperm(length(Ii))[1:new_nRBF]];
    #end
    	# println("length of I",size(Ii));
		global new_RBF_location = append!(new_RBF_location,Ii);

		if(locateRBFwithGrads)
			# println("Added rbfs at locations:",Ii);
			X1 = Ii[:,1]; X2 = Ii[:,2]; X3 = Ii[:,3];
		else
			println("ii:",Ii)
			println(ind2subv(n_tup,Ii))
			indices = getindex.(ind2subv(n_tup,Ii),[1 2 3 ]);
			println(indices)

			I1 = indices[:,1];
			I2 = indices[:,2];
			I3 = indices[:,3];

			X1 = I1*mesh.h[1] .- 0.5*mesh.h[1] .+ mesh.domain[1];
			X2 = I2*mesh.h[2] .- 0.5*mesh.h[2] .+ mesh.domain[3];
			X3 = I3*mesh.h[3] .- 0.5*mesh.h[3] .+ mesh.domain[5];
		end
		# println("X1:",X1);println("X2:",X2);println("X3:",X3);

		n_mc_old = numParamOfRBF*nRBF;
		n_mc_new = numParamOfRBF*nRBF + size(Ii,1)*numParamOfRBF;
		mc_new = zeros(length(mc) + size(Ii,1)*numParamOfRBF);
		if isRBF10
			mref = pInv.mref[:,1];
		else
			mref = pInv.mref;
		end
		mref_new = zeros(Float32,length(mc) + size(Ii,1)*numParamOfRBF);

		mc_new[1:n_mc_old] = mc[1:n_mc_old];

		mc_new[(n_mc_new+1):end] = mc[(n_mc_old+1):end];
		mref_new[(n_mc_new+1):end] = mref[(n_mc_old+1):end];

		mc_new[(n_mc_old+1):numParamOfRBF:n_mc_new] .= 0.0; #ParamLevelSet.centerHeavySide - ParamLevelSet.deltaHeavySide;
		mc_new[(n_mc_old+numParamOfRBF-2):numParamOfRBF:n_mc_new] = X1;
		mc_new[(n_mc_old+numParamOfRBF-1):numParamOfRBF:n_mc_new] = X2;
		mc_new[(n_mc_old+numParamOfRBF):numParamOfRBF:n_mc_new] = X3;

		boundsLow = zeros(Float32,length(mc_new)).+a;
		boundsHigh = zeros(Float32,length(mc_new)).+b;
		boundsLow[1:n_mc_old] = pInv.boundsLow[1:n_mc_old];
		boundsHigh[1:n_mc_old] = pInv.boundsHigh[1:n_mc_old];

		eeps =  1/(0.01171875*2.5)
		if isRBF10
			mc_new[(n_mc_old+2):10:n_mc_new] .= eeps #10.0;
			mc_new[(n_mc_old+5):10:n_mc_new] .= eeps#10.0;
			mc_new[(n_mc_old+7):10:n_mc_new] .= eeps#10.0;
			boundsLow[(n_mc_old+2):10:n_mc_new] .= 0.05;
			boundsLow[(n_mc_old+5):10:n_mc_new] .= 0.05;
			boundsLow[(n_mc_old+7):10:n_mc_new] .= 0.05;
		else
			mc_new[(n_mc_old+2):5:n_mc_new] .= 6.0;
			boundsLow[(n_mc_old+2):5:n_mc_new] .= 0.05;
		end
		# boundsLow[end-9:end] .= -0.01;
		# boundsHigh[end-9:end] .= 0.01;
		pInv.boundsLow = boundsLow;
		pInv.boundsHigh = boundsHigh;
		mref_new[1:n_mc_new] = mc_new[1:n_mc_new];
		mref_new[end-5*npc+1:end].=0.0;
		global mc = mc_new;
		global nRBF += new_nRBF;

		#II = speye(Float32,length(mc_new));
		len = Int64(length(mc_new));
		IIs = sparse(1.0I,len,len);
		pInv.regularizer = (m, mref, M)->TikhonovReg(m,mref,M,IIs);

		if isRBF10
			spd_reg = (m,mref,M) -> RBF_SPD_regularization(m,mref,nRBF);
			pInv.regularizer  = [pInv.regularizer;spd_reg];
			pInv.mref 			  = [mref_new zeros(length(mref_new))];
			pInv.alpha .*= 0.8;
		else
			pInv.mref = mref_new;
			pInv.alpha *= 0.8;
		end

		mc,Dc,flag,his = projGN(mc,pInv,pMisRFs,solveGN=projGNexplicit);
		#mc[end-5*npc+1:end] .= 0.0;
		myDump(mc,0,iterNum,pInv,0,methodName,0,noiseAnglesDeg,noiseTrans,invertVis);

		if iterNum==outerIter
			pInv.alpha *= 0.1;
			pInv.maxIter = 200;
			mc, = projGN(mc,pInv,pMisRFs,solveGN=projGNexplicit);
			myDump(mc,0,iterNum+1,pInv,0,methodName,0,noiseAnglesDeg,noiseTrans,invertVis);
		end

	end

		meshForPlotting = getRegularMesh(Mesh.domain,Mesh.n*1);
		ntup = tuple((meshForPlotting.n)...);
		if isSimple
			modfunForPlotting = (m)->MeshFreeParamLevelSetModelFunc(meshForPlotting,m,Xc = P,computeJacobian = 0,
				sigma = getDefaultHeavySide() ,bf = 1, numParamOfRBF = numParamOfRBF);
		else
			modfunForPlotting = (m)->MeshFreeParamLevelSetModelFunc(meshForPlotting,m[1:(numParamOfRBF*nRBF)],Xc = P,computeJacobian = 0,
				sigma = getDefaultHeavySide() ,bf = 1, numParamOfRBF = numParamOfRBF);
		end
		fullMc = modfunForPlotting(mc)[1];
		writedlm(string("Fandisk200Direct_Volume_final_",ntup,"_from_",tuple((Mesh.n)...),".dat"),convert(Array{Float16},fullMc));
		#m = convert(Array{Float16},m);
		#writedlm(string("Fandisk200Direct_Volume_orig.dat"),m);
		#writedlm(string("Fandisk200Direct_RBF_centers.dat"),mCenters);

end

max_stars_repo_name: BGUCompSci/ShapeReconstructionPaLS

max_issues_repo_name: BGUCompSci/ShapeReconstructionPaLS

max_forks_repo_name: BGUCompSci/ShapeReconstructionPaLS