text:
/*
	Copyright (C) 2003-2013 by David White <davewx7@gmail.com>

    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 2 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program.  If not, see <http://www.gnu.org/licenses/>.
*/
#pragma once
#if defined(USE_SHADERS)

#include <map>
#include <vector>
#include <boost/intrusive_ptr.hpp>
#include <boost/shared_ptr.hpp>

#include "asserts.hpp"
#include "entity_fwd.hpp"
#include "formula_callable.hpp"
#include "formula_callable_definition.hpp"
#include "formula_variable_storage.hpp"

namespace gles2
{

class shader
{
public:
	static void set_runtime_error(const std::string& msg);
	static std::string get_and_clear_runtime_error();

	shader() : type_(0), shader_(0)
	{}
	explicit shader(GLenum type, const std::string& name, const variant& code);
	GLuint get() const { return shader_; }
	std::string name() const { return name_; }
	std::string code() const { return code_; }
protected:
	bool compile(const std::string& code);
private:
	GLenum type_;
	GLuint shader_;
	std::string name_;
	std::string code_;
};

struct actives
{
	// Name of variable.
	std::string name;
	// type of the uniform/attribute variable
	GLenum type;
	// If an array type, this is the maximum number of array elements used
	// in the program. Value is 1 if type is not an array type.
	GLsizei num_elements;
	// Location of the active uniform/attribute
	GLint location;
	// Last value
	variant last_value;
};

class program;
typedef boost::intrusive_ptr<program> program_ptr;
typedef boost::intrusive_ptr<const program> const_program_ptr;

typedef std::map<std::string,actives>::iterator actives_map_iterator;
typedef std::map<std::string,actives>::const_iterator const_actives_map_iterator;

class program : public game_logic::formula_callable
{
public:
	program();
	explicit program(const std::string& name, const shader& vs, const shader& fs);
	virtual ~program()
	{}
	void init(const std::string& name, const shader& vs, const shader& fs);
	GLuint get() const { return object_; }
	GLuint get_attribute(const std::string& attr) const;
	GLint get_uniform(const std::string& attr) const;
	std::string name() const { return name_; }
	game_logic::formula_ptr create_formula(const variant& v);
	bool execute_command(const variant& var);
	actives_map_iterator get_uniform_reference(const std::string& name);
	actives_map_iterator get_attribute_reference(const std::string& name);
	void set_uniform(const actives_map_iterator& it, const variant& value);
	void set_uniform(const actives_map_iterator& it, const GLsizei count, const GLfloat* fv);
	void set_uniform_or_defer(actives_map_iterator& it, const variant& value);
	void set_uniform_or_defer(const std::string& key, const variant& value);
	variant get_uniform_value(const std::string& key) const;
	void set_attributes(const std::string& key, const variant& value);
	void set_attributes(const actives_map_iterator& it, const variant& value);
	variant get_attributes_value(const std::string& key) const;
	game_logic::formula_callable* get_environment() { return environ_; }
	void set_deferred_uniforms();
	GLint mvp_matrix_uniform() const { return u_mvp_matrix_; }
	GLint vertex_attribute() const { return vertex_location_; }
	GLint texcoord_attribute() const { return texcoord_location_; }
	GLuint get_fixed_attribute(const std::string& name) const;
	GLuint get_fixed_uniform(const std::string& name) const;

	virtual variant write();

	void vertex_attrib_array(GLint ndx, GLint size, GLenum type, GLboolean normalized, GLsizei stride, const GLvoid* ptr);

	virtual void vertex_array(GLint size, GLenum type, GLboolean normalized, GLsizei stride, const GLvoid* ptr);
	virtual void texture_array(GLint size, GLenum type, GLboolean normalized, GLsizei stride, const GLvoid* ptr);
	virtual void color_array(GLint size, GLenum type, GLboolean normalized, GLsizei stride, const GLvoid* ptr);
	virtual void set_fixed_attributes(const variant& node);
	virtual void set_fixed_uniforms(const variant& node);
	virtual void set_fixed_attributes();
	virtual void set_fixed_uniforms();

	void disable_vertex_attrib(GLint);

	static void load_shaders(const std::string& shader_data);
	static program_ptr find_program(const std::string& prog_name);
	static void add_shader(const std::string& program_name,
		const shader& v_shader,
		const shader& f_shader,
		const variant& prog,
		const variant& uniforms);
	static std::map<std::string, gles2::program_ptr>& get_shaders();
	static void clear_shaders();
	void set_known_uniforms();
	void set_sprite_area(const GLfloat* fl);
	void set_draw_area(const GLfloat* fl);
	void set_cycle(int cycle);

	const shader& vertex_shader() const { return vs_; }
	const shader& fragment_shader() const { return fs_; }
private:
	bool link();
	bool queryUniforms();
	bool queryAttributes();

	std::vector<GLint> active_attributes_;
	variant stored_attributes_;
	variant stored_uniforms_;

	DECLARE_CALLABLE(program)

	game_logic::formula_callable* environ_;
	std::string name_;
	shader vs_;
	shader fs_;
	GLuint object_;
	std::map<std::string, actives> attribs_;
	std::map<std::string, actives> uniforms_;

	//"vertex" and "texcoord" values within stored_attributes_
	std::string vertex_attribute_;
	std::string texcoord_attribute_;
	std::string color_attribute_;

	GLint vertex_location_, texcoord_location_, color_attr_location_;

	std::vector<std::map<std::string, actives>::iterator> uniforms_to_update_;

	GLint u_tex_map_;
	GLint u_mvp_matrix_;
	GLint u_sprite_area_;
	GLint u_draw_area_;
	GLint u_cycle_;
	GLint u_color_;
	GLint u_point_size_;
	GLint u_discard_;

	friend class shader_program;
};

class shader_program : public game_logic::formula_callable
{
public:
	static boost::intrusive_ptr<shader_program> get_global(const std::string& key);
	shader_program();
	shader_program(const shader_program& o);
	explicit shader_program(const variant& node, entity* obj = NULL);
	explicit shader_program(const std::string& program_name);
	virtual ~shader_program()
	{}

	bool enabled() const { return enabled_; }
	void configure(const variant& node, entity* obj = NULL);
	void init(entity* obj);
	game_logic::formula_ptr create_formula(const variant& v);
	bool execute_command(const variant& var);
	int zorder() const { return zorder_; }

	void prepare_draw();
	void refresh_for_draw();
	program_ptr shader() const;
	const std::string& name() const { return name_; }
	entity* parent() const { return parent_; }
	void set_parent(entity* obj) { parent_ = obj; }

	game_logic::formula_callable* vars() { return vars_.get(); }
	const game_logic::formula_callable* vars() const { return vars_.get(); }

	void clear();
	variant write();
protected:

	struct DrawCommand {
		DrawCommand();
		std::map<std::string, actives>::iterator target;
		variant value;
		bool increment;
	};

	class attribute_commands_callable : public game_logic::formula_callable
	{
	public:
		void set_program(program_ptr program) { program_ = program; }
		void execute_on_draw();
	private:
		virtual variant get_value(const std::string& key) const;
		virtual void set_value(const std::string& key, const variant& value);

		program_ptr program_;
		std::vector<DrawCommand> attribute_commands_;
	};

	class uniform_commands_callable : public game_logic::formula_callable
	{
	public:
		void set_program(program_ptr program) { program_ = program; }
		void execute_on_draw();
	private:
		virtual variant get_value(const std::string& key) const;
		virtual void set_value(const std::string& key, const variant& value);

		program_ptr program_;
		std::vector<DrawCommand> uniform_commands_;
	};

private:
	DECLARE_CALLABLE(shader_program)

	//void operator=(const shader_program&);

	std::string name_;
	program_ptr program_object_;

	game_logic::formula_variable_storage_ptr vars_;

	std::vector<std::string> create_commands_;
	std::vector<std::string> draw_commands_;
	std::vector<game_logic::formula_ptr> create_formulas_;
	std::vector<game_logic::formula_ptr> draw_formulas_;

	boost::intrusive_ptr<uniform_commands_callable> uniform_commands_;
	boost::intrusive_ptr<attribute_commands_callable> attribute_commands_;

	// fake zorder value
	int zorder_;

	entity* parent_;

	bool enabled_;
};

typedef boost::intrusive_ptr<shader_program> shader_program_ptr;
typedef boost::intrusive_ptr<const shader_program> const_shader_program_ptr;

}

GLenum get_blend_mode(variant v);

#endif

max_stars_repo_name: sweetkristas/anura

max_issues_repo_name: sweetkristas/anura

max_forks_repo_name: sweetkristas/anura


text:
#include <stdio.h>
#include <string>
#include <opencv2/opencv.hpp>
#include <ros/ros.h>
#include <cv_bridge/cv_bridge.h>
#include <sensor_msgs/Image.h>
#include <sensor_msgs/CameraInfo.h>
#include <sensor_msgs/image_encodings.h>
#include <sensor_msgs/distortion_models.h>
#include <image_transport/image_transport.h>
#include <boost/optional.hpp>
#include <boost/property_tree/ptree.hpp>
#include <boost/property_tree/ini_parser.hpp>
#include <camera_info_manager/camera_info_manager.h>

#define WIDTH_ID 3
#define HEIGHT_ID 4
#define FPS_ID 5

namespace arti {

class StereoCamera
{

public:
	/**
	 * @brief      { stereo camera driver }
	 *
	 * @param[in]  resolution  The resolution
	 * @param[in]  frame_rate  The frame rate
	 */
	 StereoCamera(int device_id, int resolution, double frame_rate):frame_rate_(30.0) {

		camera_ = new cv::VideoCapture(device_id);
		cv::Mat raw;
		cv::Mat left_image;
		cv::Mat right_image;
		setResolution(resolution);
		// // this function doesn't work very well in current Opencv 2.4, so, just use ROS to control frame rate.
		// setFrameRate(frame_rate);

		std::cout << "Stereo Camera Set Resolution: " << camera_->get(WIDTH_ID) << "x" << camera_->get(HEIGHT_ID) << std::endl;
		// std::cout << "Stereo Camera Set Frame Rate: " << camera_->get(FPS_ID) << std::endl;
	}

	~StereoCamera() {
		// std::cout << "Destroy the pointer" << std::endl;
		delete camera_;
	}

	/**
	 * @brief      Sets the resolution.
	 *
	 * @param[in]  type  The type
	 */
	void setResolution(int type) {

		if (type == 0) { width_ = 4416; height_ = 1242;} // 2K
		if (type == 1) { width_ = 3840; height_ = 1080;} // FHD
		if (type == 2) { width_ = 2560; height_ = 720;}  // HD
		if (type == 3) { width_ = 1344; height_ = 376;}  // VGA

		camera_->set(WIDTH_ID, width_);
		camera_->set(HEIGHT_ID, height_);

		// make sure that the number set are right from the hardware
		width_ = camera_->get(WIDTH_ID);
		height_ = camera_->get(HEIGHT_ID);

	}

	/**
	 * @brief      Sets the frame rate.
	 *
	 * @param[in]  frame_rate  The frame rate
	 */
	void setFrameRate(double frame_rate) {
		camera_->set(FPS_ID, frame_rate);
		frame_rate_ = camera_->get(FPS_ID);
	}

	/**
	 * @brief      Gets the images.
	 *
	 * @param      left_image   The left image
	 * @param      right_image  The right image
	 *
	 * @return     The images.
	 */
	bool getImages(cv::Mat& left_image, cv::Mat& right_image) {
		cv::Mat raw;
		if (camera_->grab()) {
			camera_->retrieve(raw);
			cv::Rect left_rect(0, 0, width_ / 2, height_);
			cv::Rect right_rect(width_ / 2, 0, width_ / 2, height_);
			left_image = raw(left_rect);
			right_image = raw(right_rect);
			cv::waitKey(10);
			return true;
		} else {
			return false;
		}
	}

private:
	cv::VideoCapture* camera_;
	int width_;
	int height_;
	double frame_rate_;
	bool cv_three_;
};

/**
 * @brief       the camera ros warpper class
 */
class ZedCameraROS {
public:

	/**
	 * @brief      { function_description }
	 *
	 * @param[in]  resolution  The resolution
	 * @param[in]  frame_rate  The frame rate
	 */
	ZedCameraROS() {
		ros::NodeHandle nh;
		ros::NodeHandle private_nh("~");
		// get ros param
		private_nh.param("resolution", resolution_, 1);
		private_nh.param("frame_rate", frame_rate_, 30.0);
		private_nh.param("config_file_location", config_file_location_, std::string("~/SN1267.conf"));
		private_nh.param("left_frame_id", left_frame_id_, std::string("left_camera"));
		private_nh.param("right_frame_id", right_frame_id_, std::string("right_camera"));
		private_nh.param("show_image", show_image_, false);
		private_nh.param("load_zed_config", load_zed_config_, true);
		private_nh.param("device_id", device_id_, 0);

		correctFramerate(resolution_,frame_rate_);

		ROS_INFO("Try to initialize the camera");
		StereoCamera zed(device_id_, resolution_, frame_rate_);
		ROS_INFO("Initialized the camera");

		// setup publisher stuff
		image_transport::ImageTransport it(nh);
		image_transport::Publisher left_image_pub = it.advertise("left/image_raw", 1);
		image_transport::Publisher right_image_pub = it.advertise("right/image_raw", 1);

		ros::Publisher left_cam_info_pub = nh.advertise<sensor_msgs::CameraInfo>("left/camera_info", 1);
		ros::Publisher right_cam_info_pub = nh.advertise<sensor_msgs::CameraInfo>("right/camera_info", 1);

		sensor_msgs::CameraInfo left_info, right_info;

		ROS_INFO("Try load camera calibration files");
		if (load_zed_config_) {
			ROS_INFO("Loading from zed calibration files");
			// get camera info from zed
			try {
				getZedCameraInfo(config_file_location_, resolution_, left_info, right_info);
			}
			catch (std::runtime_error& e) {
				ROS_INFO("Can't load camera info");
				ROS_ERROR("%s", e.what());
				throw e;
			}
		} else {
			ROS_INFO("Loading from ROS calibration files");
			// get config from the left, right.yaml in config
			camera_info_manager::CameraInfoManager info_manager(nh);
			info_manager.setCameraName("zed/left");
			info_manager.loadCameraInfo( "package://zed_cpu_ros/config/left.yaml");
			left_info = info_manager.getCameraInfo();

			info_manager.setCameraName("zed/right");
			info_manager.loadCameraInfo( "package://zed_cpu_ros/config/right.yaml");
			right_info = info_manager.getCameraInfo();

			left_info.header.frame_id = left_frame_id_;
			right_info.header.frame_id = right_frame_id_;
		}

		// std::cout << left_info << std::endl;
		// std::cout << right_info << std::endl;

		ROS_INFO("Got camera calibration files");
		// loop to publish images;
		cv::Mat left_image, right_image;
		ros::Rate r(frame_rate_);

		while (nh.ok()) {
			ros::Time now = ros::Time::now();
			if (!zed.getImages(left_image, right_image)) {
				ROS_INFO_ONCE("Can't find camera");
			} else {
				ROS_INFO_ONCE("Success, found camera");
			}
			if (show_image_) {
				cv::imshow("left", left_image);
				cv::imshow("right", right_image);
			}

			if (left_image_pub.getNumSubscribers() > 0) {
				publishImage(left_image, left_image_pub, "left_frame", now);
			}
			if (right_image_pub.getNumSubscribers() > 0) {
				publishImage(right_image, right_image_pub, "right_frame", now);
			}
			if (left_cam_info_pub.getNumSubscribers() > 0) {
				publishCamInfo(left_cam_info_pub, left_info, now);
			}
			if (right_cam_info_pub.getNumSubscribers() > 0) {
				publishCamInfo(right_cam_info_pub, right_info, now);
			}
			r.sleep();
			// since the frame rate was set inside the camera, no need to do a ros sleep
		}
	}

	/**
	 * @brief      Gets the camera information From Zed config.
	 *
	 * @param[in]  config_file         The configuration file
	 * @param[in]  resolution          The resolution
	 * @param[in]  left_cam_info_msg   The left camera information message
	 * @param[in]  right_cam_info_msg  The right camera information message
	 */
	void getZedCameraInfo(std::string config_file, int resolution, sensor_msgs::CameraInfo& left_info, sensor_msgs::CameraInfo& right_info) {
		boost::property_tree::ptree pt;
		boost::property_tree::ini_parser::read_ini(config_file, pt);
		std::string left_str = "LEFT_CAM_";
		std::string right_str = "RIGHT_CAM_";
		std::string reso_str = "";

		switch (resolution) {
			case 0: reso_str = "2K"; break;
			case 1: reso_str = "FHD"; break;
			case 2: reso_str = "HD"; break;
			case 3: reso_str = "VGA"; break;
		}
		// left value
		double l_cx = pt.get<double>(left_str + reso_str + ".cx");
		double l_cy = pt.get<double>(left_str + reso_str + ".cy");
		double l_fx = pt.get<double>(left_str + reso_str + ".fx");
		double l_fy = pt.get<double>(left_str + reso_str + ".fy");
		double l_k1 = pt.get<double>(left_str + reso_str + ".k1");
		double l_k2 = pt.get<double>(left_str + reso_str + ".k2");
		// right value
		double r_cx = pt.get<double>(right_str + reso_str + ".cx");
		double r_cy = pt.get<double>(right_str + reso_str + ".cy");
		double r_fx = pt.get<double>(right_str + reso_str + ".fx");
		double r_fy = pt.get<double>(right_str + reso_str + ".fy");
		double r_k1 = pt.get<double>(right_str + reso_str + ".k1");
		double r_k2 = pt.get<double>(right_str + reso_str + ".k2");

		// get baseline and convert mm to m
		boost::optional<double> baselineCheck;
		double baseline = 0.0;
		// some config files have "Baseline" instead of "BaseLine", check accordingly...
		if (baselineCheck = pt.get_optional<double>("STEREO.BaseLine")) {
			baseline = pt.get<double>("STEREO.BaseLine") * 0.001;
		}
		else if (baselineCheck = pt.get_optional<double>("STEREO.Baseline")) {
			baseline = pt.get<double>("STEREO.Baseline") * 0.001;
		}
		else
			throw std::runtime_error("baseline parameter not found");

		// get Rx and Rz
		double rx = pt.get<double>("STEREO.RX_"+reso_str);
		double rz = pt.get<double>("STEREO.RZ_"+reso_str);
		double ry = pt.get<double>("STEREO.CV_"+reso_str);

		// assume zeros, maybe not right
		double p1 = 0, p2 = 0, k3 = 0;

		left_info.distortion_model = sensor_msgs::distortion_models::PLUMB_BOB;
		right_info.distortion_model = sensor_msgs::distortion_models::PLUMB_BOB;

		// distortion parameters
		// For "plumb_bob", the 5 parameters are: (k1, k2, t1, t2, k3).
		left_info.D.resize(5);
		left_info.D[0] = l_k1;
		left_info.D[1] = l_k2;
		left_info.D[2] = k3;
		left_info.D[3] = p1;
		left_info.D[4] = p2;

		right_info.D.resize(5);
		right_info.D[0] = r_k1;
		right_info.D[1] = r_k2;
		right_info.D[2] = k3;
		right_info.D[3] = p1;
		right_info.D[4] = p2;

		// Intrinsic camera matrix
		// 	[fx  0 cx]
		// K =  [ 0 fy cy]
		//	[ 0  0  1]
		left_info.K.fill(0.0);
		left_info.K[0] = l_fx;
		left_info.K[2] = l_cx;
		left_info.K[4] = l_fy;
		left_info.K[5] = l_cy;
		left_info.K[8] = 1.0;

		right_info.K.fill(0.0);
		right_info.K[0] = r_fx;
		right_info.K[2] = r_cx;
		right_info.K[4] = r_fy;
		right_info.K[5] = r_cy;
		right_info.K[8] = 1.0;

		// rectification matrix
		// Rl = R_rect, R_r = R * R_rect
		// since R is identity, Rl = Rr;
		left_info.R.fill(0.0);
		right_info.R.fill(0.0);
		cv::Mat rvec = (cv::Mat_<double>(3, 1) << rx, ry, rz);
		cv::Mat rmat(3, 3, CV_64F);
		cv::Rodrigues(rvec, rmat);
		int id = 0;
		cv::MatIterator_<double> it, end;
		for (it = rmat.begin<double>(); it != rmat.end<double>(); ++it, id++){
			left_info.R[id] = *it;
			right_info.R[id] = *it;
		}

		// Projection/camera matrix
		//     [fx'  0  cx' Tx]
		// P = [ 0  fy' cy' Ty]
		//     [ 0   0   1   0]
		left_info.P.fill(0.0);
		left_info.P[0] = l_fx;
		left_info.P[2] = l_cx;
		left_info.P[5] = l_fy;
		left_info.P[6] = l_cy;
		left_info.P[10] = 1.0;

		right_info.P.fill(0.0);
		right_info.P[0] = r_fx;
		right_info.P[2] = r_cx;
		right_info.P[3] = (-1 * l_fx * baseline);
		right_info.P[5] = r_fy;
		right_info.P[6] = r_cy;
		right_info.P[10] = 1.0;

		left_info.width = right_info.width = width_;
		left_info.height = right_info.height = height_;

		left_info.header.frame_id = left_frame_id_;
		right_info.header.frame_id = right_frame_id_;
	}

	/**
	 * @brief      { publish camera info }
	 *
	 * @param[in]  pub_cam_info  The pub camera information
	 * @param[in]  cam_info_msg  The camera information message
	 * @param[in]  now           The now
	 */
	void publishCamInfo(const ros::Publisher& pub_cam_info, sensor_msgs::CameraInfo& cam_info_msg, ros::Time now) {
		cam_info_msg.header.stamp = now;
		pub_cam_info.publish(cam_info_msg);
	}

	/**
	 * @brief      { publish image }
	 *
	 * @param[in]  img           The image
	 * @param      img_pub       The image pub
	 * @param[in]  img_frame_id  The image frame identifier
	 * @param[in]  t             { parameter_description }
	 */
	void publishImage(cv::Mat img, image_transport::Publisher &img_pub, std::string img_frame_id, ros::Time t) {
		cv_bridge::CvImage cv_image;
		cv_image.image = img;
		cv_image.encoding = sensor_msgs::image_encodings::BGR8;
		cv_image.header.frame_id = img_frame_id;
		cv_image.header.stamp = t;
		img_pub.publish(cv_image.toImageMsg());
	}
	/**
	 * @brief      Correct frame rate according to resolution
	 *
	 * @param[in]  resolution          The resolution
	 * @param      frame_rate   			 The camera frame rate
	 */
	 void correctFramerate(int resolution, double &frame_rate){
		 double max_frame_rate;
		 std::string reso_str = "";
		 switch (resolution) {
 			case 0: max_frame_rate = 15; reso_str = "2K"; break;
 			case 1: max_frame_rate = 30; reso_str = "FHD"; break;
 			case 2: max_frame_rate = 60; reso_str = "HD"; break;
 			case 3: max_frame_rate = 100; reso_str = "VGA"; break;
 		}
		if(frame_rate > max_frame_rate)
			ROS_WARN("frame_rate(%fHz) to high for resolution(%s), downgraded to %fHz",frame_rate,reso_str.c_str(),max_frame_rate);
			frame_rate = max_frame_rate;
	 }

private:
	int device_id_, resolution_;
	double frame_rate_;
	bool show_image_, load_zed_config_;
	double width_, height_;
	std::string left_frame_id_, right_frame_id_;
	std::string config_file_location_;
};

}

int main(int argc, char **argv) {
    try {
        ros::init(argc, argv, "zed_camera");
        arti::ZedCameraROS zed_ros;
        return EXIT_SUCCESS;
    }
    catch(std::runtime_error& e) {
        ros::shutdown();
        return EXIT_FAILURE;
    }
}

ext: cpp


text:
/*
 *  Copyright Matus Chochlik.
 *  Distributed under the Boost Software License, Version 1.0.
 *  See accompanying file LICENSE_1_0.txt or copy at
 *   http://www.boost.org/LICENSE_1_0.txt
 */

#define BOOST_TEST_DYN_LINK
#include "../random.hpp"
#include <eagine/diagnostic.hpp>
#include <eagine/interop/valgrind.hpp>
#include <eagine/maybe_unused.hpp>

#ifdef __clang__
EAGINE_DIAG_PUSH()
EAGINE_DIAG_OFF(disabled-macro-expansion)
#endif

#include <boost/test/tools/floating_point_comparison.hpp>
#include <boost/test/unit_test.hpp>

namespace eagine {
//------------------------------------------------------------------------------
template <typename T>
static inline T test_repeats(T min, T max) noexcept {
    return running_on_valgrind() ? min : max;
}
//------------------------------------------------------------------------------
} // namespace eagine

ext: inl

max_issues_repo_name: matus-chochlik/eagine-core

max_forks_repo_name: matus-chochlik/eagine-core


text:
#define _FILE_OFFSET_BITS 64
#undef _FORTIFY_SOURCE

#include <reader.h>
#include <images.h>
#include <timing.h>

#include <string.h>
#include <stdlib.h>
#include <stdio.h>
#ifdef _WIN32
  #include <windows.h>
#else
  #include <sys/mman.h>
#endif
#include <math.h>

#include <set>
#include <map>
#include <vector>

#include <boost/bind.hpp>
#include <boost/function.hpp>

using namespace std;

struct circuit_map {
  int *data;
  int nl, sx, sy;

  int p(int l, int x, int y) const {
    return data[nl*(x+y*sx)+l];
  }

  void s(int l, int x, int y, int v) {
    data[nl*(x+y*sx)+l] = v;
  }

  circuit_map(const char *fname, int nl, int sx, int sy, bool create);
  ~circuit_map();
};

circuit_map::circuit_map(const char *fname, int _nl, int _sx, int _sy, bool create)
{
  nl = _nl;
  sx = _sx;
  sy = _sy;
  unsigned char *map_adr;
  if(create) {
    create_file_rw(fname, map_adr, nl*4*(int64_t)sx*sy);
    memset(map_adr, 0xff, (int64_t)sx*sy*4*nl);

  } else {
    int64_t size;
    map_file_ro(fname, map_adr, size, false);
  }

  data = (int *)map_adr;
}

circuit_map::~circuit_map()
{
  #ifdef _WIN32
    UnmapViewOfFile(data);
  #else
    munmap(data, long(nl)*sx*sy*4);
  #endif
}

enum {
  ACTIVE,
  POLY,
  METAL,
  BURIED,
  TRANSISTOR,
  CAPACITOR,
  DISABLED
};

const char *type_names = "apmbtc*";

struct circuit_info {
  int type;
  int surface;
  int x0, y0, x1, y1;
  int net, netp, metal;
  std::set<int> neighbors;
};

struct via_info {
  int metal;
  int active_poly;

  via_info(int _metal, int _active_poly) {
    metal = _metal;
    active_poly = _active_poly;
  }
};

struct net_info {
  set<int> circuits;
};

struct trans_info {
  int circ, t1, t2, gate, x, y;
  double strength;
};

struct via_map {
  std::map<int, std::vector<int> > metal_to_ap, ap_to_metal;
};

struct fill_coord {
  int x, y;
  fill_coord(int _x, int _y) { x = _x; y = _y; }
};

struct metal_link_info {
  int x1, y1, x2, y2;
  int circuit_id;
};

void fill_1(std::vector<fill_coord> &fc, int x, int y, int sx, int sy, int color, boost::function<void(int, int)> set_pixel, boost::function<int(int, int)> read_color, boost::function<bool(int, int)> test_pixel)
{
  if(test_pixel(x, y))
    return;
  set_pixel(x, y);
  int x0 = x-1;
  while(x0 >= 0 && !test_pixel(x0, y) && read_color(x0, y) == color) {
    set_pixel(x0, y);
    x0--;
  }
  x0++;
  int x1 = x+1;
  while(x1 < sx && !test_pixel(x1, y) && read_color(x1, y) == color) {
    set_pixel(x1, y);
    x1++;
  }
  x1--;
  for(int xx=x0; xx<=x1; xx++) {
    if(y > 0   && !test_pixel(xx, y-1) && read_color(xx, y-1) == color)
      fc.push_back(fill_coord(xx, y-1));
    if(y < sy-1 && !test_pixel(xx, y+1) && read_color(xx, y+1) == color)
      fc.push_back(fill_coord(xx, y+1));
  }
}

void fill(int x, int y, int sx, int sy, int color, boost::function<void(int, int)> set_pixel, boost::function<int(int, int)> read_color, boost::function<bool(int, int)> test_pixel)
{
  std::vector<fill_coord> fc;
  fc.push_back(fill_coord(x, y));
  while(!fc.empty()) {
    int xx = fc.front().x;
    int yy = fc.front().y;
    fc.erase(fc.begin());
    fill_1(fc, xx, yy, sx, sy, color, set_pixel, read_color, test_pixel);
  }
}

void color_active_poly_set(int x, int y, int *pmap, bool *is_buried, bool *is_caps)
{
  int v = pmap[x+y*11];
  if(v & 2)
    *is_buried = true;
  if(v & 4)
    *is_caps = true;
  pmap[x+y*11] = v | 8;
}

int color_active_poly_read(int x, int y, int *pmap)
{
  return pmap[x+y*11] & 1;
}

bool color_active_poly_test(int x, int y, int *pmap)
{
  return pmap[x+y*11] & 8;
}

int color_active_poly(int x, int y, const pbm *active, const pbm *poly, const pbm *buried, const pbm *caps)
{
  bool ca = active->p(x, y);
  bool cp = poly->p(x, y);
  if(ca && cp) {
    for(int yy = y-5; yy <= y+5; yy++)
      for(int xx = x-5; xx <= x+5; xx++)
	if(buried->p(xx, yy) || (caps && caps->p(xx, yy)))
	  goto check_buried_caps;
    return TRANSISTOR;
  check_buried_caps:
    {
      int pmap[11*11];
      for(int yy = 0; yy <= 10; yy++)
	for(int xx = 0; xx <= 10; xx++) {
	  int x1 = x-5+xx;
	  int y1 = y-5+yy;
	  pmap[yy*11+xx] = active->p(x1, y1) && poly->p(x1, y1) ? buried->p(x1, y1) ? 3 : caps && caps->p(x1, y1) ? 5 : 1 : 0;
	}
      bool is_buried = false;
      bool is_caps = false;
      fill(5, 5, 11, 11, 1, boost::bind(color_active_poly_set, _1, _2, pmap, &is_buried, & is_caps), boost::bind(color_active_poly_read, _1, _2, pmap), boost::bind(color_active_poly_test, _1, _2, pmap));
      return is_caps ? CAPACITOR : is_buried ? BURIED : TRANSISTOR;
    }
  }

  return ca ? ACTIVE : cp ? POLY : -1;
}

int color_active_gates(int x, int y, const pbm *active, const pbm *gates, const pbm *caps)
{
  bool ca = active->p(x, y);
  bool cg = gates->p(x, y);
  bool cc = caps && caps->p(x, y);
  assert(!(cg && cc));
  if(cc)
    return CAPACITOR;
  if(ca)
    return ACTIVE;
  if(cg)
    return TRANSISTOR;
  return -1;
}

int color_metal(int x, int y, const pbm *metal)
{
  return metal->p(x, y) ? METAL : -1;
}

void circuit_set(int x, int y, int cid, circuit_info *ci, circuit_map *dest, int l)
{
  if(ci->x0 > x)
    ci->x0 = x;
  if(ci->y0 > y)
    ci->y0 = y;
  if(ci->x1 < x)
    ci->x1 = x;
  if(ci->y1 < y)
    ci->y1 = y;
  ci->surface++;
  dest->s(l, x, y, cid);
}

bool circuit_test(int x, int y, int cid, const circuit_map *dest, int l)
{
  return dest->p(l, x, y) == cid;
}

void build_circuits(time_info &tinfo, boost::function<int(int, int)> color, std::vector<circuit_info> &circuits, circuit_map *dest, int l)
{
  for(int y=0; y<dest->sy; y++) {
    tinfo.tick(y, dest->sy);
    for(int x=0; x<dest->sx; x++) {
      if(dest->p(l, x, y) != -1)
	continue;
      int c = color(x, y);
      if(c == -1)
	continue;
      int cid = circuits.size();
      circuits.resize(cid+1);
      circuit_info *ci = &circuits[cid];
      ci->type = c;
      ci->surface = 0;
      ci->x0 = ci->x1 = x;
      ci->y0 = ci->y1 = y;
      ci->net = -1;
      ci->netp = -1;
      ci->metal = -1;
      fill(x, y, dest->sx, dest->sy, c, boost::bind(circuit_set, _1, _2, cid, ci, dest, l), color, boost::bind(circuit_test, _1, _2, cid, dest, l));
    }
  }
}

static void build_neighbors_add(std::vector<circuit_info> &circuits, int cm, int x, int y, int cn)
{
  if(cn == -1 || cn == cm)
    return;
  int cmt = circuits[cm].type;
  int cnt = circuits[cn].type;
  if((cmt == ACTIVE && cnt != POLY) || (cmt == POLY && cnt != ACTIVE) || (cmt != POLY && cmt != ACTIVE))
    circuits[cm].neighbors.insert(cn);
}

void build_neighbors(time_info &tinfo, std::vector<circuit_info> &circuits, const circuit_map &cmap)
{
  for(int y=1; y<cmap.sy-1; y++) {
    tinfo.tick(y-1, cmap.sy-2);
    for(int x=1; x<cmap.sx-1; x++) {
      int cm = cmap.p(0, x, y);
      if(cm == -1)
	continue;
      build_neighbors_add(circuits, cm, x-1, y, cmap.p(0, x-1, y));
      build_neighbors_add(circuits, cm, x+1, y, cmap.p(0, x+1, y));
      build_neighbors_add(circuits, cm, x, y-1, cmap.p(0, x, y-1));
      build_neighbors_add(circuits, cm, x, y+1, cmap.p(0, x, y+1));
    }
  }
}

void build_transistor_groups(std::vector<set<int> > &groups, int &terminaux, int &gates, std::vector<bool> &is_terminal, std::vector<bool> &is_gate, const std::vector<circuit_info> &circuits, int id)
{
  const circuit_info &ci = circuits[id];
  terminaux = 0;
  gates = 0;
  for(int step=0; step<2; step++) {
    for(set<int>::const_iterator j = ci.neighbors.begin(); j != ci.neighbors.end(); j++) {
      const circuit_info &ci1 = circuits[*j];
      if(step ^ (ci1.type != BURIED && ci1.type != CAPACITOR))
	continue;
      std::vector<set<int> >::iterator l;
      for(l = groups.begin(); l != groups.end(); l++) {
	for(set<int>::const_iterator k = l->begin(); k != l->end(); k++)
	  if(ci1.neighbors.find(*k) != ci1.neighbors.end()) {
	    const circuit_info &ci2 = circuits[*k];
	    if(!((ci1.type == ACTIVE && ci2.type == POLY) || (ci1.type == POLY && ci2.type == ACTIVE)))
	      goto found;
	  }
      }

    found:
      if(l == groups.end()) {
	groups.push_back(set<int>());
	groups.back().insert(*j);
	if(ci1.type == ACTIVE || ci1.type == BURIED || ci1.type == CAPACITOR) {
	  terminaux++;
	  is_terminal.push_back(true);
	} else
	  is_terminal.push_back(false);

	if(ci1.type == POLY || ci1.type == BURIED || ci1.type == CAPACITOR) {
	  gates++;
	  is_gate.push_back(true);
	} else
	  is_gate.push_back(false);

      } else
	l->insert(*j);
    }
  }
}

void clean_and_remap(time_info &tinfo, std::vector<circuit_info> &circuits, circuit_map &cmap)
{
  bool has_error = false;

  int nc = circuits.size();
  int *remap_active = new int[nc];
  int *remap_poly   = new int[nc];

  for(unsigned int i=0; i != circuits.size(); i++) {
    tinfo.tick(i, circuits.size());
    circuit_info &ci = circuits[i];
    switch(ci.type) {
    case ACTIVE:
      remap_active[i] = i;
      remap_poly[i] = -1;
      break;
    case POLY:
      remap_active[i] = -1;
      remap_poly[i] = i;
      break;
    case BURIED:
      remap_active[i] = i;
      remap_poly[i] = i;
      break;
    case CAPACITOR:
      remap_active[i] = i;
      remap_poly[i] = i;
      break;
    case METAL:
      remap_active[i] = -1;
      remap_poly[i] = -1;
      break;
    case TRANSISTOR: {
      std::vector<set<int> > groups;
      std::vector<bool> is_gate, is_terminal;
      int terminaux, gates;
      build_transistor_groups(groups, terminaux, gates, is_terminal, is_gate, circuits, i);
      if(terminaux < 2) { //  || (ci.x0 - ci.x1 <= 2 && ci.x0 - ci.x1 >= -2) || (ci.y0 - ci.y1 <= 2 && ci.y0 - ci.y1 >= -2)) {
	if(terminaux == 0 || gates == 0) {
	  fprintf(stderr, "P/A superposition zone (%d, %d)-(%d, %d) has no active %s\n", ci.x0, cmap.sy-1-ci.y1, ci.x1, cmap.sy-1-ci.y0, terminaux ? "gate" : "terminal");
	  for(set<int>::const_iterator j = ci.neighbors.begin(); j != ci.neighbors.end(); j++) {
	    const circuit_info &ci1 = circuits[*j];
	    fprintf(stderr, "  - %c%d (%d, %d)-(%d, %d)\n", type_names[ci1.type], *j, ci1.x0, cmap.sy-1-ci1.y1, ci1.x1, cmap.sy-1-ci1.y0);
	  }
	  has_error = true;
	} else {
	  int pmap = -1, amap = -1;
	  for(set<int>::const_iterator j = ci.neighbors.begin(); j != ci.neighbors.end(); j++) {
	    int t = circuits[*j].type;
	    switch(t) {
	    case ACTIVE: amap = *j; break;
	    case POLY:   pmap = *j; break;
	    case BURIED: if(amap == -1) amap = *j; if(pmap == -1) pmap = *j; break;
	    }
	  }
	  if(amap == -1 || pmap == -1) {
	    fprintf(stderr, "Bad a/pmap on clean-and-remap (%d, %d)-(%d, %d)\n", ci.x0, cmap.sy-1-ci.y1, ci.x1, cmap.sy-1-ci.y0);
	    has_error = true;
	    break;
	  }
	  for(set<int>::const_iterator j = ci.neighbors.begin(); j != ci.neighbors.end(); j++) {
	    circuits[amap].neighbors.insert(*j);
	    circuits[pmap].neighbors.insert(*j);
	  }

	  ci.type = DISABLED;
	  ci.neighbors.clear();
	  remap_active[i] = amap;
	  remap_poly[i] = pmap;
	}
      } else {
	remap_active[i] = i;
	remap_poly[i] = i;
      }
      break;
    }
    }
  }

  if(has_error)
    exit(1);

  for(int y=0; y<cmap.sy; y++) {
    tinfo.tick(y, cmap.sy);
    for(int x=0; x<cmap.sx; x++) {
      int v = cmap.p(0, x, y);
      if(v != -1) {
	cmap.s(0, x, y, remap_active[v]);
	cmap.s(1, x, y, remap_poly[v]);
      }
    }
  }

  for(unsigned int i=0; i != circuits.size(); i++) {
    circuit_info &ci = circuits[i];
  retry:
    for(set<int>::iterator j = ci.neighbors.begin(); j != ci.neighbors.end(); j++)
      if(remap_active[*j] != -1 && remap_active[*j] != *j) {
	int rm = ci.type == ACTIVE ? remap_active[*j] : remap_poly[*j];
	ci.neighbors.erase(j);
	ci.neighbors.insert(rm);
	goto retry;
      }
  }
}

void compress_ids(time_info &tinfo, std::vector<circuit_info> &circuit_infos, std::vector<metal_link_info> &virtual_poly_id, circuit_map &cmap)
{
  std::vector<int> remap;
  remap.resize(circuit_infos.size());
  unsigned int cid = 0;
  for(unsigned int i=0; i != circuit_infos.size(); i++) {
    remap[i] = cid;
    if(circuit_infos[i].type != DISABLED) {
      if(cid != i)
	circuit_infos[cid] = circuit_infos[i];
      cid++;
    }
  }
  circuit_infos.resize(cid);

  for(int y=0; y<cmap.sy; y++) {
    tinfo.tick(y, cmap.sy);
    for(int x=0; x<cmap.sx; x++)
      for(int l=0; l<3; l++) {
	int v = cmap.p(l, x, y);
	if(v != -1)
	  cmap.s(l, x, y, remap[v]);
      }
  }

  for(std::vector<metal_link_info>::iterator i = virtual_poly_id.begin(); i != virtual_poly_id.end(); i++)
    i->circuit_id = remap[i->circuit_id];

  for(unsigned int i=0; i != circuit_infos.size(); i++) {
    circuit_info &ci = circuit_infos[i];
    set<int> re;
    for(set<int>::const_iterator i = ci.neighbors.begin(); i != ci.neighbors.end(); i++)
      re.insert(remap[*i]);
    ci.neighbors = re;
  }
}

void map_vias_set(int x, int y, via_info *via, const std::vector<circuit_info> *circuit_infos, pbm *used, circuit_map *cmap)
{
  used->s(x, y, true);
  int na = cmap->p(0, x, y);
  int np = cmap->nl >= 3 ? cmap->p(1, x, y) : -1;
  int nm = cmap->p(cmap->nl >= 3 ? 2 : 1, x, y);
  if(nm != -1) {
    if(via->metal == -1)
      via->metal = nm;
    else if(via->metal != nm)
      fprintf(stderr, "via at (%d, %d) touches multiple metal tracks\n", x, cmap->sy-1-y);
  }
  if(na != -1 && np != -1 && na != np)
    fprintf(stderr, "via at (%d, %d) touches split active/poly\n", x, cmap->sy-1-y);
  else if(na != -1 || np != -1) {
    int nn = na == -1 ? np : na;
    if((*circuit_infos)[nn].type == TRANSISTOR)
      fprintf(stderr, "via at (%d, %d) touches a transistor\n", x, cmap->sy-1-y);
    else if((*circuit_infos)[nn].type == CAPACITOR)
      fprintf(stderr, "via at (%d, %d) touches a capacitor\n", x, cmap->sy-1-y);
    else if(via->active_poly == -1)
      via->active_poly = nn;
    else if(via->active_poly != nn)
      fprintf(stderr, "via at (%d, %d) touches multiple poly/layer zones\n", x, cmap->sy-1-y);
  }
}

void map_vias(time_info &tinfo, std::vector<via_info> &via_infos, via_map &via_maps, const std::vector<circuit_info> &circuit_infos, const pbm *vias, circuit_map &cmap)
{
  pbm used(cmap.sx, cmap.sy);
  for(int y=0; y<cmap.sy; y++) {
    tinfo.tick(y, cmap.sy);
    for(int x=0; x<cmap.sx; x++) {
      if(vias->p(x, y) && !used.p(x, y)) {
	via_info via(-1, -1);
	fill(x, y, cmap.sx, cmap.sy, 1, boost::bind(map_vias_set, _1, _2, &via, &circuit_infos, &used, &cmap), boost::bind(&pbm::p, vias, _1, _2), boost::bind(&pbm::p, &used, _1, _2));
	if(via.metal == -1)
	  fprintf(stderr, "via at (%d, %d) does not touch the metal\n", x, cmap.sy-1-y);
	if(via.active_poly == -1)
	  fprintf(stderr, "via at (%d, %d) does not touch poly or active\n", x, cmap.sy-1-y);

	int yy = cmap.sy-1-y;
	//	if(!(x >= 6000 && x < 6750 && yy >= 5000 && yy < 6500))
	//	  continue;
	if(0 && x > 6500 && x < 7000 && yy >= 6500 && yy < 7000)
	  continue;
	if(via.metal != -1 && via.active_poly != -1)
	  via_infos.push_back(via);
      }
    }
  }

  for(unsigned int i=0; i != via_infos.size(); i++) {
    via_maps.metal_to_ap[via_infos[i].metal].push_back(via_infos[i].active_poly);
    via_maps.ap_to_metal[via_infos[i].active_poly].push_back(via_infos[i].metal);
  }
}

void via_list_add(std::vector<int> &stack, const map<int, std::vector<int> > vmap, int id, const std::vector<circuit_info> &circuit_infos)
{
  map<int, std::vector<int> >::const_iterator vi = vmap.find(id);
  if(vi == vmap.end())
    return;
  for(std::vector<int>::const_iterator i = vi->second.begin(); i != vi->second.end(); i++)
    if(circuit_infos[*i].net == -1)
      stack.push_back(*i);
}

void build_nets(time_info &tinfo, std::vector<net_info> &net_infos, std::vector<circuit_info> &circuit_infos, const via_map &via_maps, bool has_poly, int sx, int sy)
{
  for(unsigned int i=0; i != circuit_infos.size(); i++) {
    //    tinfo.tick(i, circuit_infos.size());
    if(circuit_infos[i].type != DISABLED && circuit_infos[i].net == -1) {
      if(circuit_infos[i].type == CAPACITOR || (circuit_infos[i].type == TRANSISTOR && !has_poly))
	continue;
      int nid = net_infos.size();
      net_infos.resize(nid+1);
      set<int> &circuits = net_infos[nid].circuits;
      std::vector<int> stack;
      assert(circuit_infos[i].type != CAPACITOR);
      stack.push_back(i + (circuit_infos[i].type == POLY ? 1000000 : 0));
      do {
	int cid = stack.front();
	int pcid = cid;
	stack.erase(stack.begin());
	int start_stack_size = stack.size();
	int is_poly = cid / 1000000;
	cid = cid % 1000000;
	circuit_info &ci = circuit_infos[cid];
	int curnet = is_poly && ci.type == CAPACITOR ? ci.netp : ci.net;
	if(curnet != -1) {
	  if(curnet != nid) {
	    fprintf(stderr, "Network creation failure on %c%d (%d, %d)-(%d, %d), nets %d and %d want to link\n", type_names[ci.type], cid, ci.x0, sy-1-ci.y1, ci.x1, sy-1-ci.y0, nid, curnet);
	    fprintf(stderr, "  pcid=%d\n", pcid);
	    fprintf(stderr, "  net %d (%d/%d):\n", curnet, ci.net, ci.netp);
	    for(set<int>::const_iterator j = net_infos[curnet].circuits.begin(); j != net_infos[curnet].circuits.end(); j++) {
	      const circuit_info &ci1 = circuit_infos[*j];
	      fprintf(stderr, "  - %c%d (%d/%d) (%d, %d)-(%d, %d)", type_names[ci1.type], *j, ci.net, ci.netp, ci1.x0, sy-1-ci1.y1, ci1.x1, sy-1-ci1.y0);
	      for(set<int>::const_iterator j = ci1.neighbors.begin(); j != ci1.neighbors.end(); j++) {
		const circuit_info &ci2 = circuit_infos[*j];
		fprintf(stderr, " %c%d", type_names[ci2.type], *j);
	      }
	      fprintf(stderr, "\n");
	    }
	    fprintf(stderr, "  net %d:\n", nid);
	    for(set<int>::const_iterator j = net_infos[nid].circuits.begin(); j != net_infos[nid].circuits.end(); j++) {
	      const circuit_info &ci1 = circuit_infos[*j];
	      fprintf(stderr, "  - %c%d (%d, %d)-(%d, %d)", type_names[ci1.type], *j, ci1.x0, sy-1-ci1.y1, ci1.x1, sy-1-ci1.y0);
	      for(set<int>::const_iterator j = ci1.neighbors.begin(); j != ci1.neighbors.end(); j++) {
		const circuit_info &ci2 = circuit_infos[*j];
		fprintf(stderr, " %c%d", type_names[ci2.type], *j);
	      }
	      fprintf(stderr, "\n");
	    }
	    exit(1);
	  }
	  continue;
	}

	circuits.insert(cid + (ci.type == CAPACITOR && is_poly ? 1000000 : 0));
	if(ci.type == CAPACITOR && is_poly)
	  ci.netp = nid;
	else
	  ci.net = nid;
	switch(ci.type) {
	case ACTIVE:
	  via_list_add(stack, via_maps.ap_to_metal, cid, circuit_infos);
	  for(set<int>::const_iterator j = ci.neighbors.begin(); j != ci.neighbors.end(); j++)
	    if(circuit_infos[*j].type == BURIED || circuit_infos[*j].type == CAPACITOR)
	      stack.push_back(*j);
	  break;

	case POLY:
	  via_list_add(stack, via_maps.ap_to_metal, cid, circuit_infos);
	  for(set<int>::const_iterator j = ci.neighbors.begin(); j != ci.neighbors.end(); j++) {
	    int id = *j;
	    if(circuit_infos[id].net == -1 && (circuit_infos[id].type == BURIED || circuit_infos[id].type == TRANSISTOR))
	      stack.push_back(id + 1000000);
	    if(circuit_infos[id].netp == -1 && circuit_infos[id].type == CAPACITOR)
	      stack.push_back(id + 1000000);
	  }
	  break;

	case METAL:
	  via_list_add(stack, via_maps.metal_to_ap, cid, circuit_infos);
	  break;

	case BURIED:
	  via_list_add(stack, via_maps.ap_to_metal, cid, circuit_infos);
	  for(set<int>::const_iterator j = ci.neighbors.begin(); j != ci.neighbors.end(); j++) {
	    if(circuit_infos[*j].type == ACTIVE || circuit_infos[*j].type == CAPACITOR)
	      stack.push_back(*j);
	    if(circuit_infos[*j].type == POLY || circuit_infos[*j].type == TRANSISTOR)
	      stack.push_back(*j + 1000000);
	    if(circuit_infos[*j].type == CAPACITOR)
	      stack.push_back(*j + 1000000);
	  }
	  break;

	case TRANSISTOR:
	  for(set<int>::const_iterator j = ci.neighbors.begin(); j != ci.neighbors.end(); j++) {
	    if(circuit_infos[*j].type == BURIED || circuit_infos[*j].type == POLY || circuit_infos[*j].type == CAPACITOR)
	      stack.push_back(*j + 1000000);
	  }
	  break;

	case CAPACITOR:
	  if(is_poly) {
	    for(set<int>::const_iterator j = ci.neighbors.begin(); j != ci.neighbors.end(); j++)
	      if(circuit_infos[*j].type == BURIED || circuit_infos[*j].type == POLY || circuit_infos[*j].type == TRANSISTOR)
		stack.push_back(*j + 1000000);
	  } else {
	    for(set<int>::const_iterator j = ci.neighbors.begin(); j != ci.neighbors.end(); j++)
	      if(circuit_infos[*j].type == BURIED || circuit_infos[*j].type == ACTIVE)
		stack.push_back(*j);
	  }
	  break;

	default:
	  abort();
	}

	if(false && !nid) {
	  std::vector<int>::const_iterator i = stack.begin();
	  for(int ii=0; ii<start_stack_size; ii++)
	    i++;
	  printf("%d:", cid);
	  for(;i != stack.end();i++)
	    printf(" %d", *i);
	  printf("\n");
	}
#if 0
	fprintf(stderr, "  + %c%d %d (%d, %d)-(%d, %d)", type_names[ci.type], cid, is_poly, ci.x0, cmap.sy-1-ci.y1, ci.x1, cmap.sy-1-ci.y0);
	fprintf(stderr, "    ");
	for(std::vector<int>::const_iterator j = stack.begin(); j != stack.end(); j++)
	  fprintf(stderr, " %d", *j);
	fprintf(stderr, "\n");
#endif
      } while(!stack.empty());
    }
  }
}

static double dist_min(double v1, double v2)
{
  if(v2 >= 0 && (v1 < 0 || v2 < v1))
    return v2;
  else
    return v1;
}

double build_transistors_compute_force_dir(const set<int> &src, const set<int> &dst, int gate, const circuit_info &gci, const circuit_map &cmap, int &middle_x, int &middle_y)
{
  int sx = gci.x1-gci.x0+1 + 4;
  int sy = gci.y1-gci.y0+1 + 4;
  double *grid = new double[sx*sy];
  std::vector<pair<int, int> > dests;
  int x0 = gci.x0 - 2;
  int y0 = gci.y0 - 2;
  double *g = grid;
  for(int y=0; y<sy; y++)
    for(int x=0; x<sx; x++) {
      int cet = cmap.p(0, x+x0, y+y0);
      double v = -100000;
      if(cet == gate)
	v = 100000;
      else if(src.find(cet) != src.end())
	v = 0;
      *g++ = v;
      if(x != 0 && x != sx-1 && y != 0 && y != sy-1 && dst.find(cet) != dst.end())
	dests.push_back(pair<int, int>(x, y));
    }

  int minx=1, miny=1, maxx=sx-2, maxy=sy-2;
  while(minx<=maxx) {
    int minx1 = minx;
    int miny1 = miny;
    int maxx1 = maxx;
    int maxy1 = maxy;
    minx = sx;
    miny = sy;
    maxx = -1;
    maxy = -1;
    for(int y=miny1; y<=maxy1; y++)
      for(int x=minx1; x<=maxx1; x++) {
	g = grid + y*sx + x;
	if(*g <= 0)
	  continue;
	double v = -100000;
	v = dist_min(v, g[-1-sx] + M_SQRT2);
	v = dist_min(v, g[  -sx] + 1);
	v = dist_min(v, g[ 1-sx] + M_SQRT2);
	v = dist_min(v, g[-1   ] + 1);
	v = dist_min(v, g[ 1   ] + 1);
	v = dist_min(v, g[-1+sx] + M_SQRT2);
	v = dist_min(v, g[   sx] + 1);
	v = dist_min(v, g[ 1+sx] + M_SQRT2);
	if(v < *g) {
	  *g = v;
	  if(x-1 < minx)
	    minx = x-1;
	  if(x+1 > maxx)
	    maxx = x+1;
	  if(y-1 < miny)
	    miny = y-1;
	  if(y+1 > maxy)
	    maxy = y+1;
	}
      }
  }
  double f = 0;
  double midx = 0, midy = 0;
  for(std::vector<pair<int, int> >::const_iterator i = dests.begin(); i != dests.end(); i++) {
    g = grid + i->second*sx + i->first;
    double v = -100000;
    v = dist_min(v, g[-1-sx] + M_SQRT2);
    v = dist_min(v, g[  -sx] + 1);
    v = dist_min(v, g[ 1-sx] + M_SQRT2);
    v = dist_min(v, g[-1   ] + 1);
    v = dist_min(v, g[ 1   ] + 1);
    v = dist_min(v, g[-1+sx] + M_SQRT2);
    v = dist_min(v, g[   sx] + 1);
    v = dist_min(v, g[ 1+sx] + M_SQRT2);
    if(v > 0) {
      double ff = 1/v;
      f += ff;
      midx += i->first * ff;
      midy += i->second * ff;
    }
  }
  middle_x = x0 + int(midx / f + 0.5);
  middle_y = y0 + int(midy / f + 0.5);
  delete[] grid;
  return f;
}

double build_transistors_compute_force(const set<int> &t1, const set<int> &t2, int gate, const circuit_info &gci, const circuit_map &cmap, int &middle_x, int &middle_y)
{
  double f1, f2;
  int x1, x2, y1, y2;
  f1 = build_transistors_compute_force_dir(t1, t2, gate, gci, cmap, x1, y1);
  f2 = build_transistors_compute_force_dir(t2, t1, gate, gci, cmap, x2, y2);
  middle_x = (x1+x2)/2;
  middle_y = (y1+y2)/2;

  return 0.5*(f1 + f2);
}

void build_transistors(time_info &tinfo, std::vector<trans_info> &trans_infos, const std::vector<net_info> &net_infos, const std::vector<circuit_info> &circuit_infos, const circuit_map &cmap)
{
  for(unsigned int i=0; i != circuit_infos.size(); i++) {
    tinfo.tick(i, circuit_infos.size());
    if(circuit_infos[i].type == TRANSISTOR) {
      std::vector<set<int> > groups;
      std::vector<bool> is_gate, is_terminal;
      int terminaux, gates;
      build_transistor_groups(groups, terminaux, gates, is_terminal, is_gate, circuit_infos, i);

      int ng = groups.size();
      double *forces = new double[ng*ng];
      int *midx = new int[ng*ng];
      int *midy = new int[ng*ng];
      int j=0;
      for(std::vector<set<int> >::const_iterator ji = groups.begin(); ji != groups.end(); ji++, j++) {
	if(!is_terminal[j])
	  continue;
	int k = j+1;
	std::vector<set<int> >::const_iterator ki = ji;
	ki++;
	for(;ki != groups.end();ki++, k++) {
	  if(!is_terminal[k])
	    continue;
	  int idx = j*ng+k;
	  forces[idx] = build_transistors_compute_force(*ji, *ki, i, circuit_infos[i], cmap, midx[idx], midy[idx]);
	}
      }

      std::vector<int> nets;
      std::vector<bool> linked;
      linked.resize(ng);
      for(j=0; j<ng; j++)
	linked[j] = false;

      for(std::vector<set<int> >::const_iterator ji = groups.begin(); ji != groups.end(); ji++)
	nets.push_back(circuit_infos[*ji->begin()].net);

      for(j=0; j<ng; j++) {
	if(linked[j] || !is_terminal[j])
	  continue;
	int best_k = -1;
	double best_f = 0;
	int best_idx = 0;
	for(int k=0; k<ng; k++) {
	  if(j == k || !is_terminal[k])
	    continue;
	  int idx = j < k ? j*ng+k : k*ng+j;
	  double f = forces[idx];
	  if(best_k == -1 || f > best_f) {
	    best_k = k;
	    best_f = f;
	    best_idx = idx;
	  }
	}
	linked[j] = true;
	linked[best_k] = true;
	trans_info ti;
	ti.circ = i;
	ti.t1 = nets[j];
	ti.t2 = nets[best_k];
	ti.gate = circuit_infos[i].net;
	ti.strength = best_f;
	ti.x = midx[best_idx];
	ti.y = midy[best_idx];
	trans_infos.push_back(ti);
      }
      delete[] forces;
      delete[] midx;
      delete[] midy;
    }
  }
}

void build_transistors_metal_gate(time_info &tinfo, std::vector<trans_info> &trans_infos, const std::vector<net_info> &net_infos, const std::vector<circuit_info> &circuit_infos, const circuit_map &cmap)
{
  for(unsigned int i=0; i != circuit_infos.size(); i++) {
    tinfo.tick(i, circuit_infos.size());
    if(circuit_infos[i].type == TRANSISTOR) {
      const circuit_info &ci = circuit_infos[i];
      if(ci.neighbors.size() != 2) {
	fprintf(stderr, "Gate at (%d, %d)-(%d, %d) has %d neighbors.\n",
		ci.x0, cmap.sy-1-ci.y0, ci.x1, cmap.sy-1-ci.y1,
		int(ci.neighbors.size()));
	continue;
      }

      trans_info ti;
      ti.circ = i;

      set<int>::const_iterator ni = ci.neighbors.begin();
      set<int> t1, t2;
      t1.insert(*ni);
      ti.t1 = circuit_infos[*ni++].net;
      t2.insert(*ni);
      ti.t2 = circuit_infos[*ni++].net;
      ti.gate = ci.net;
      ti.strength = build_transistors_compute_force(t1, t2, i, ci, cmap, ti.x, ti.y);
      trans_infos.push_back(ti);
    }
  }
}

void circuit_stats(const std::vector<circuit_info> &circuit_infos)
{
  int t[7];
  t[0] = t[1] = t[2] = t[3] = t[4] = t[5] = t[6] = 0;
  for(unsigned int i=0; i != circuit_infos.size(); i++)
    t[circuit_infos[i].type]++;

  fprintf(stderr, "-> %d active, %d poly, %d metal, %d buried, %d gates, %d capacitors, %d disabled, %d total\n",
	  t[0], t[1], t[2], t[3], t[4], t[5], t[6], t[0]+t[1]+t[2]+t[3]+t[4]+t[5]+t[6]);
}

void dump(FILE *out, int sx, int sy, int nl, const std::vector<circuit_info> &circuit_infos)
{
  fprintf(out, "%d %d %d\n", sx, sy, nl);
  fprintf(out, "%d circuits\n", int(circuit_infos.size()));
  for(unsigned int i=0; i != circuit_infos.size(); i++) {
    const circuit_info &ci = circuit_infos[i];
    fprintf(out, "%6d %c %5d %5d %5d %5d %5d %5d %d", i, type_names[ci.type], ci.net, ci.netp, ci.x0, sy-1-ci.y1, ci.x1, sy-1-ci.y0, ci.surface);
    for(set<int>::const_iterator j = ci.neighbors.begin(); j != ci.neighbors.end(); j++)
      fprintf(out, " %c%d", type_names[circuit_infos[*j].type], *j);
    fprintf(out, "\n");
  }
}

void dump(FILE *out, int sx, int sy, int nl, const std::vector<net_info> &net_infos)
{
  fprintf(out, "%d nets\n", int(net_infos.size()));
  for(unsigned int i=0; i != net_infos.size(); i++) {
    const net_info &ni = net_infos[i];
    fprintf(out, "%5d", i);
    for(set<int>::const_iterator j = ni.circuits.begin(); j != ni.circuits.end(); j++)
      fprintf(out, " %d", *j);
    fprintf(out, "\n");
  }
}

void dump(FILE *out, int sx, int sy, int nl, const std::vector<trans_info> &trans_infos)
{
  fprintf(out, "%d transistors\n", int(trans_infos.size()));
  for(unsigned int i=0; i != trans_infos.size(); i++) {
    const trans_info &ti = trans_infos[i];
    fprintf(out, "%5d %6d %5d %5d %5d %5d %5d %g\n", i, ti.circ, ti.x, sy-1 - ti.y, ti.t1, ti.gate, ti.t2, ti.strength);
  }
}

void dump(const char *fname, int sx, int sy, int nl, const std::vector<trans_info> &trans_infos, const std::vector<net_info> &net_infos, const std::vector<circuit_info> &circuit_infos)
{
  FILE *out = fopen(fname, "w");
  dump(out, sx, sy, nl, circuit_infos);
  dump(out, sx, sy, nl, net_infos);
  dump(out, sx, sy, nl, trans_infos);
  fclose(out);
}

// Add a virtual circuit between the two vcc domains
void add_virtual_circuit(std::vector<circuit_info> &circuit_infos, metal_link_info &ml, int sy, bool has_poly)
{
  int id = circuit_infos.size();
  circuit_infos.resize(id+1);
  circuit_info &ci = circuit_infos[id];
  ci.type = has_poly ? POLY :  ACTIVE;
  ci.surface = abs(ml.x2-ml.x1)+abs(ml.y2-ml.y1)+1;
  ci.x0 = ml.x1;
  ci.x1 = ml.x2;
  ci.y1 = sy-1-ml.y2;
  ci.y0 = sy-1-ml.y1;
  ci.net = -1;
  ml.circuit_id = id;
}

// Add a pair of virtual vias with the virtual circuit
void add_virtual_vias(std::vector<via_info> &via_infos, via_map &via_maps, metal_link_info &ml, const circuit_map &cmap)
{
  int vcc1 = cmap.p(cmap.nl >= 3 ? 2 : 1, ml.x1, cmap.sy-1-ml.y1);
  int vcc2 = cmap.p(cmap.nl >= 3 ? 2 : 1, ml.x2, cmap.sy-1-ml.y2);
  via_infos.push_back(via_info(vcc1, ml.circuit_id));
  via_infos.push_back(via_info(vcc2, ml.circuit_id));
  via_maps.metal_to_ap[vcc1].push_back(ml.circuit_id);
  via_maps.metal_to_ap[vcc2].push_back(ml.circuit_id);
  via_maps.ap_to_metal[ml.circuit_id].push_back(vcc1);
  via_maps.ap_to_metal[ml.circuit_id].push_back(vcc2);
}

// In a metal gates circuit, lookup the metal net corresponding to gates and caps
void lookup_gates_and_caps(time_info &tinfo, std::vector<circuit_info> &circuit_infos, const circuit_map &cmap)
{
  for(int y=0; y<cmap.sy; y++) {
    tinfo.tick(y, cmap.sy);
    for(int x=0; x<cmap.sx; x++) {
      int ca = cmap.p(0, x, y);
      if(ca == -1)
	continue;
      circuit_info &ci = circuit_infos[ca];
      if(ci.type == TRANSISTOR || ci.type == CAPACITOR) {
	int cm = cmap.p(1, x, y);
	if(cm == -1) {
	  fprintf(stderr, "%s does not touch metal at (%d, %d)\n", ci.type == TRANSISTOR ? "Gate" : "Capacitor", x, cmap.sy-1-y);
	  exit(1);
	}
	int nm = circuit_infos[cm].net;
	if(ci.type == TRANSISTOR) {
	  if(ci.net == -1)
	    ci.net = nm;
	  else if(ci.net != nm) {
	    fprintf(stderr, "Transistor touches multiple metal at (%d, %d)\n", x, cmap.sy-1-y);
	    exit(1);
	  }
	} else {
	  if(ci.netp == -1)
	    ci.netp = nm;
	  else if(ci.netp != nm) {
	    fprintf(stderr, "Capacitor touches multiple metal at (%d, %d)\n", x, cmap.sy-1-y);
	    exit(1);
	  }
	}
      }
    }
  }
}

pbm *active = NULL;
pbm *gates = NULL;
pbm *buried = NULL;
pbm *metal = NULL;
pbm *poly = NULL;
pbm *vias = NULL;
pbm *caps = NULL;
std::vector<metal_link_info> metal_links;

const char *map_name = NULL;
const char *list_name = NULL;
int sx;
int sy;

void nmos_poly_single_metal()
{
  std::vector<circuit_info> circuit_infos;
  std::vector<via_info> via_infos;
  via_map via_maps;
  std::vector<net_info> net_infos;
  std::vector<trans_info> trans_infos;

  circuit_map cmap(map_name, 3, sx, sy, true);

  time_info tinfo;
  tinfo.start("build circuits active/poly");
  build_circuits(tinfo, boost::bind(color_active_poly, _1, _2, active, poly, buried, caps), circuit_infos, &cmap, 0);
  for(std::vector<metal_link_info>::iterator i = metal_links.begin(); i != metal_links.end(); i++)
    add_virtual_circuit(circuit_infos, *i, sy, true);
  circuit_stats(circuit_infos);
  tinfo.start("build circuits metal");
  build_circuits(tinfo, boost::bind(color_metal, _1, _2, metal), circuit_infos, &cmap, 2);
  circuit_stats(circuit_infos);
  tinfo.start("build neighbors");
  build_neighbors(tinfo, circuit_infos, cmap);
  tinfo.start("clean and remap");
  clean_and_remap(tinfo, circuit_infos, cmap);
  circuit_stats(circuit_infos);
  tinfo.start("compressing ids");
  compress_ids(tinfo, circuit_infos, metal_links, cmap);
  circuit_stats(circuit_infos);
  tinfo.start("mapping vias");
  map_vias(tinfo, via_infos, via_maps, circuit_infos, vias, cmap);
  for(std::vector<metal_link_info>::iterator i = metal_links.begin(); i != metal_links.end(); i++)
    add_virtual_vias(via_infos, via_maps, *i, cmap);
  fprintf(stderr, "  -> %d vias mapped\n", int(via_infos.size()));
  tinfo.start("building nets");
  build_nets(tinfo, net_infos, circuit_infos, via_maps, true, sx, sy);
  fprintf(stderr, "  -> %d nets built\n", int(net_infos.size()));
  tinfo.start("building transistors");
  build_transistors(tinfo, trans_infos, net_infos, circuit_infos, cmap);
  fprintf(stderr, "  -> %d transistors built\n", int(trans_infos.size()));
  fprintf(stderr, "Dumping info...\n");
  dump(list_name, sx, sy, 3, trans_infos, net_infos, circuit_infos);
}

void nmos_metal_gate()
{
  std::vector<circuit_info> circuit_infos;
  std::vector<via_info> via_infos;
  via_map via_maps;
  std::vector<net_info> net_infos;
  std::vector<trans_info> trans_infos;

  circuit_map cmap(map_name, 2, sx, sy, true);
  time_info tinfo;
  tinfo.start("build circuits active/gates");
  build_circuits(tinfo, boost::bind(color_active_gates, _1, _2, active, gates, caps), circuit_infos, &cmap, 0);
  for(std::vector<metal_link_info>::iterator i = metal_links.begin(); i != metal_links.end(); i++)
    add_virtual_circuit(circuit_infos, *i, sy, false);
  circuit_stats(circuit_infos);
  tinfo.start("build circuits metal");
  build_circuits(tinfo, boost::bind(color_metal, _1, _2, metal), circuit_infos, &cmap, 1);
  circuit_stats(circuit_infos);
  tinfo.start("build neighbors");
  build_neighbors(tinfo, circuit_infos, cmap);
  tinfo.start("mapping vias");
  map_vias(tinfo, via_infos, via_maps, circuit_infos, vias, cmap);
  for(std::vector<metal_link_info>::iterator i = metal_links.begin(); i != metal_links.end(); i++)
    add_virtual_vias(via_infos, via_maps, *i, cmap);
  fprintf(stderr, "  -> %d vias mapped\n", int(via_infos.size()));
  tinfo.start("building nets");
  build_nets(tinfo, net_infos, circuit_infos, via_maps, false, sx, sy);
  fprintf(stderr, "  -> %d nets built\n", int(net_infos.size()));
  tinfo.start("lookup gates and caps");
  lookup_gates_and_caps(tinfo, circuit_infos, cmap);
  tinfo.start("building transistors");
  build_transistors_metal_gate(tinfo, trans_infos, net_infos, circuit_infos, cmap);
  fprintf(stderr, "  -> %d transistors built\n", int(trans_infos.size()));
  fprintf(stderr, "Dumping info...\n");
  dump(list_name, sx, sy, 2, trans_infos, net_infos, circuit_infos);
}

int main(int argc, char **argv)
{
  if(argc != 2) {
    fprintf(stderr, "Usage:\n%s config.txt\n", argv[0]);
    exit(1);
  }

  reader rd(argv[1]);
  map_name = rd.gw();
  list_name = rd.gw();
  sx = rd.gi();
  sy = rd.gi();
  rd.nl();

  void (*method)() = NULL;

  while(!rd.eof()) {
    char buf[4096];
    string keyw = rd.gw();

    if(keyw == "nmos-poly-single-metal") {
      method = nmos_poly_single_metal;
      rd.nl();

    } else if(keyw == "nmos-metal-gate") {
      method = nmos_metal_gate;
      rd.nl();

    } else if(keyw == "active") {
      sprintf(buf, "%s.pbm", rd.gwnl());
      rd.nl();
      active = new pbm(buf);

    } else if(keyw == "poly") {
      sprintf(buf, "%s.pbm", rd.gwnl());
      rd.nl();
      poly = new pbm(buf);

    } else if(keyw == "metal") {
      sprintf(buf, "%s.pbm", rd.gwnl());
      rd.nl();
      metal = new pbm(buf);

    } else if(keyw == "buried") {
      sprintf(buf, "%s.pbm", rd.gwnl());
      rd.nl();
      buried = new pbm(buf);

    } else if(keyw == "vias") {
      sprintf(buf, "%s.pbm", rd.gwnl());
      rd.nl();
      vias = new pbm(buf);

    } else if(keyw == "caps") {
      sprintf(buf, "%s.pbm", rd.gwnl());
      rd.nl();
      caps = new pbm(buf);

    } else if(keyw == "gates") {
      sprintf(buf, "%s.pbm", rd.gwnl());
      rd.nl();
      gates = new pbm(buf);

    } else if(keyw == "metal-link") {
      metal_link_info ml;
      ml.x1 = rd.gi();
      ml.y1 = rd.gi();
      ml.x2 = rd.gi();
      ml.y2 = rd.gi();
      ml.circuit_id = -1;
      rd.nl();
      metal_links.push_back(ml);

    } else {
      fprintf(stderr, "Unhandled keyword %s\n", keyw.c_str());
      exit(1);
    }
  }

  if(method)
    method();
  else
    fprintf(stderr, "Method missing, nothing to do.\n");

  return 0;
}

ext: cc


text:
#pragma once

#include <polyfem/ElementBases.hpp>
#include <polyfem/AssemblerUtils.hpp>
#include <polyfem/LocalBoundary.hpp>
#include <polyfem/Mesh.hpp>

#include <polyfem/Common.hpp>

#include <vector>
#include <Eigen/Dense>
#include <memory>

namespace polyfem
{
	class Problem
	{
	public:
		Problem(const std::string &name);
		virtual void init(const Mesh &mesh) {}

		virtual void rhs(const AssemblerUtils &assembler, const std::string &formulation, const Eigen::MatrixXd &pts, const double t, Eigen::MatrixXd &val) const = 0;
		virtual bool is_rhs_zero() const = 0;

		virtual void bc(const Mesh &mesh, const Eigen::MatrixXi &global_ids, const Eigen::MatrixXd &uv, const Eigen::MatrixXd &pts, const double t, Eigen::MatrixXd &val) const = 0;
		virtual void velocity_bc(const Mesh &mesh, const Eigen::MatrixXi &global_ids, const Eigen::MatrixXd &uv, const Eigen::MatrixXd &pts, const double t, Eigen::MatrixXd &val) const {}
		virtual void acceleration_bc(const Mesh &mesh, const Eigen::MatrixXi &global_ids, const Eigen::MatrixXd &uv, const Eigen::MatrixXd &pts, const double t, Eigen::MatrixXd &val) const {}

		virtual void neumann_bc(const Mesh &mesh, const Eigen::MatrixXi &global_ids, const Eigen::MatrixXd &uv, const Eigen::MatrixXd &pts, const Eigen::MatrixXd &normals, const double t, Eigen::MatrixXd &val) const {}
		virtual void neumann_velocity_bc(const Mesh &mesh, const Eigen::MatrixXi &global_ids, const Eigen::MatrixXd &uv, const Eigen::MatrixXd &pts, const Eigen::MatrixXd &normals, const double t, Eigen::MatrixXd &val) const {}
		virtual void neumann_acceleration_bc(const Mesh &mesh, const Eigen::MatrixXi &global_ids, const Eigen::MatrixXd &uv, const Eigen::MatrixXd &pts, const Eigen::MatrixXd &normals, const double t, Eigen::MatrixXd &val) const {}

		virtual void exact(const Eigen::MatrixXd &pts, const double t, Eigen::MatrixXd &val) const {};
		virtual void exact_grad(const Eigen::MatrixXd &pts, const double t, Eigen::MatrixXd &val) const {};

		inline const std::string &name() const { return name_; }

		virtual bool has_exact_sol() const = 0;
		virtual bool is_scalar() const = 0;
		// virtual bool is_mixed() const { return false; }
		virtual bool is_constant_in_time() const { return true; }

		virtual void clear() {}

		virtual bool is_time_dependent() const { return false; }
		virtual void initial_solution(const Mesh &mesh, const Eigen::MatrixXi &global_ids, const Eigen::MatrixXd &pts, Eigen::MatrixXd &val) const {}
		virtual void initial_velocity(const Mesh &mesh, const Eigen::MatrixXi &global_ids, const Eigen::MatrixXd &pts, Eigen::MatrixXd &val) const {}
		virtual void initial_acceleration(const Mesh &mesh, const Eigen::MatrixXi &global_ids, const Eigen::MatrixXd &pts, Eigen::MatrixXd &val) const {}
		virtual void initial_density(const Eigen::MatrixXd &pts, Eigen::MatrixXd &val) const {}

		virtual int n_incremental_load_steps(const double diag) const { return 1; }

		virtual void set_parameters(const json &params) {}

		void setup_bc(const Mesh &mesh, const std::vector<ElementBases> &bases, const std::vector<ElementBases> &pressure_bases, std::vector<LocalBoundary> &local_boundary, std::vector<int> &boundary_nodes, std::vector<LocalBoundary> &local_neumann_boundary, std::vector<int> &pressure_boundary_nodes);

		// std::vector<int> &boundary_ids() { return boundary_ids_; }
		// const std::vector<int> &boundary_ids() const { return boundary_ids_; }

		virtual bool is_dimension_dirichet(const int tag, const int dim) const { return true; }

		//here for efficiency reasons
		virtual bool all_dimensions_dirichlet() const { return true; }

		virtual ~Problem() {}

	protected:
		std::vector<int> boundary_ids_;
		std::vector<int> neumann_boundary_ids_;
		std::vector<int> pressure_boundary_ids_;
		std::vector<int> splitting_pressure_boundary_ids_;

	private:
		std::string name_;
	};

	class ProblemFactory
	{
	public:
		static const ProblemFactory &factory();

		std::shared_ptr<Problem> get_problem(const std::string &problem) const;
		inline const std::vector<std::string> &get_problem_names() const { return problem_names_; }

	private:
		ProblemFactory();
		std::map<std::string, std::function<std::shared_ptr<Problem>()>> problems_;
		std::vector<std::string> problem_names_;
	};
} // namespace polyfem

ext: hpp

max_stars_repo_name: danielepanozzo/polyfem

max_issues_repo_name: danielepanozzo/polyfem

max_forks_repo_name: danielepanozzo/polyfem


text: #include <boost/preprocessor/seq/filter.hpp>

max_stars_repo_name: miathedev/BoostForArduino

max_issues_repo_name: miathedev/BoostForArduino

max_forks_repo_name: miathedev/BoostForArduino


text:
/**
   \file gsl_simplex.hpp
   \brief a simple wrapper for the GSL_SIMPLEX method
   \author Junhua Gu
 */

#ifndef GSL_SIMPLEX_METHOD
#define GSL_SIMPLEX_METHOD
#define OPT_HEADER
#include <core/optimizer.hpp>
//#include <blitz/array.h>
#include <vector>
#include <limits>
#include <cassert>
#include <cmath>
#include <algorithm>
#include <gsl/gsl_multimin.h>
#include <iostream>

namespace opt_utilities
{

  /**
     \brief object function of the gsl simplex function
   */
  template <typename rT,typename pT>
  double gsl_func_adapter(const gsl_vector* v,void* params)
  {
    pT temp;
    temp.resize(v->size);
    for(size_t i=0;i<get_size(temp);++i)
      {
	set_element(temp,i,gsl_vector_get(v,i));
      }
    return ((func_obj<rT,pT>*)params)->eval(temp);
  }

  /**
     \brief wrapper for the gsl simplex optimization method
     \tparam return type of the object function
     \tparam param type of the object function
  */
  template <typename rT,typename pT>
  class gsl_simplex
    :public opt_method<rT,pT>
  {
  public:
    typedef pT array1d_type;
    typedef rT T;
  private:
    func_obj<rT,pT>* p_fo;
    optimizer<rT,pT>* p_optimizer;

    //typedef blitz::Array<rT,2> array2d_type;

  private:
    array1d_type start_point;
    array1d_type end_point;

  private:
    rT threshold;
  private:
    rT func(const pT& x)
    {
      assert(p_fo!=0);
      return p_fo->eval(x);
    }

    const char* do_get_type_name()const
    {
      return "gsl simplex";
    }
  public:
    gsl_simplex()
      :threshold(1e-4)
    {}

    virtual ~gsl_simplex()
    {
    };

    gsl_simplex(const gsl_simplex<rT,pT>& rhs)
      :p_fo(rhs.p_fo),p_optimizer(rhs.p_optimizer),
       start_point(rhs.start_point),
       end_point(rhs.end_point),
       threshold(rhs.threshold)
    {
    }

    gsl_simplex<rT,pT>& operator=(const gsl_simplex<rT,pT>& rhs)
    {
      threshold=rhs.threshold;
      p_fo=rhs.p_fo;
      p_optimizer=rhs.p_optimizer;
      opt_eq(start_point,rhs.start_point);
      opt_eq(end_point,rhs.end_point);
    }

    opt_method<rT,pT>* do_clone()const
    {
      return new gsl_simplex<rT,pT>(*this);
    }

    void do_set_start_point(const array1d_type& p)
    {
      start_point.resize(get_size(p));
      opt_eq(start_point,p);

    }

    array1d_type do_get_start_point()const
    {
      return start_point;
    }

    void do_set_precision(rT t)
    {
      threshold=t;
    }

    rT do_get_precision()const
    {
      return threshold;
    }

    void do_set_optimizer(optimizer<rT,pT>& o)
    {
      p_optimizer=&o;
      p_fo=p_optimizer->ptr_func_obj();
    }

    pT do_optimize()
    {
      const gsl_multimin_fminimizer_type *T =
	gsl_multimin_fminimizer_nmsimplex;
      gsl_multimin_fminimizer *s = NULL;
      gsl_vector *ss, *x;
      gsl_multimin_function minex_func;

      size_t iter = 0;
      int status;
      double size;

      /* Starting point */
      x = gsl_vector_alloc (get_size(start_point));
      //      gsl_vector_set (x, 0, 5.0);
      //gsl_vector_set (x, 1, 7.0);
      for(size_t i=0;i!=get_size(start_point);++i)
	{
	  gsl_vector_set(x,i,get_element(start_point,i));
	}

      /* Set initial step sizes to 1 */
      ss = gsl_vector_alloc (get_size(start_point));
      gsl_vector_set_all (ss, 1.0);

      //foo f;
      /* Initialize method and iterate */
      minex_func.n = get_size(start_point);
      minex_func.f = &gsl_func_adapter<double,std::vector<double> >;
      minex_func.params = (void *)p_fo;

      s = gsl_multimin_fminimizer_alloc (T, get_size(start_point));
      gsl_multimin_fminimizer_set (s, &minex_func, x, ss);

      do
	{
	  iter++;
	  status = gsl_multimin_fminimizer_iterate(s);

	  if (status)
	    {
	      break;
	    }
	  //std::cerr<<"threshold="<<threshold<<std::endl;
	  size = gsl_multimin_fminimizer_size (s);
	  status = gsl_multimin_test_size (size, threshold);

	  if (status == GSL_SUCCESS)
	    {
	      //printf ("converged to minimum at\n");
	    }

	  //printf ("%5d %10.3e %10.3ef f() = %7.3f size = %.3f\n",
	  //iter,
	  //gsl_vector_get (s->x, 0),
	  //gsl_vector_get (s->x, 1),
	  //  s->fval, size);
	}
      while (status == GSL_CONTINUE);

      /*
	foo f;
	gsl_vector_set (x, 0, 0.0);
	gsl_vector_set (x, 1, 0.0);
	cout<<"fdsa ";
	cout<<gsl_func_adapter<double,vector<double> >(x,(void*)&f)<<endl;;

      */

      end_point.resize(get_size(start_point));
      for(size_t i=0;i<get_size(start_point);++i)
	{
	  set_element(end_point,i,gsl_vector_get(s->x,i));
	}

      gsl_vector_free(x);
      gsl_vector_free(ss);
      gsl_multimin_fminimizer_free (s);

      return end_point;
    }
  };

}

#endif
//EOF


text:
#include "postgres.h"

#include <boost/lexical_cast.hpp>
#include <endian.h>
#include <stdlib.h>
#include <cstdio>
#include <cstring>
#include <map>
#include <utility>
#include <stdint.h>

#define __STDC_FORMAT_MACROS
#include <inttypes.h>

#include "../common/formatterExceptions.h"
#include "attribute.h"
#include "postgresAttribute.h"

#define __STDC_FORMAT_MACROS

Postgres::~Postgres() {
}

Postgres::Postgres(const std::string & fileName, const char * fileMode) :
    Format(fileName, fileMode) {
    setTypeAttributeMap();
}

Postgres::Postgres() :
    Format() {
}

void Postgres::setTypeAttributeMap() {
    cleanTypeAttributeMap();
    //printf("set type attribute map for Postgres\n");
    /* int32_t type */
    typeAttributeMap.insert(
        std::make_pair("int32_t",
                       new PostgresAttribute<int32_t>(fp, false)));
    typeAttributeMap.insert(
        std::make_pair("int32_t" + nullString,
                       new PostgresAttribute<int32_t>(fp, true)));
    /* int type */
    typeAttributeMap.insert(
        std::make_pair("int", new PostgresAttribute<int32_t>(fp, false)));
    typeAttributeMap.insert(
        std::make_pair("int" + nullString,
                       new PostgresAttribute<int32_t>(fp, true)));
    /* int64 type */
    typeAttributeMap.insert(
        std::make_pair("int64", new PostgresAttribute<int64_t>(fp, false)));
    typeAttributeMap.insert(
        std::make_pair("int64" + nullString,
                       new PostgresAttribute<int64_t>(fp, true)));
    /* string type */
    typeAttributeMap.insert(
        std::make_pair("string", new PostgresAttribute<char>(fp, false)));
    typeAttributeMap.insert(
        std::make_pair("string" + nullString,
                       new PostgresAttribute<char>(fp, true)));
    /* bool type */
    typeAttributeMap.insert(
        std::make_pair("bool", new PostgresAttribute<bool>(fp, false)));
    typeAttributeMap.insert(
        std::make_pair("bool" + nullString,
                       new PostgresAttribute<bool>(fp, true)));
    /* double type */
    typeAttributeMap.insert(
        std::make_pair("double", new PostgresAttribute<double>(fp, false)));
    typeAttributeMap.insert(
        std::make_pair("double" + nullString,
                       new PostgresAttribute<double>(fp, true)));
}

static const char BinarySignature[19] =
{ 'P', 'G', 'C', 'O', 'P', 'Y', '\n', '\377', '\r', '\n', '\0', '\0',
  '\0', '\0', '\0', '\0', '\0', '\0', '\0' };

bool Postgres::isTheEnd() {
    // read 2 bytes representing number of columns stored in each line
    // we have to read the colNumber and check if it is the end of the file
    // -1 represents the end of the binary data
    int16_t colNumber;
    size_t numberOfObjectsRead = fread(&colNumber, 2, 1, fp);
    if (numberOfObjectsRead != 1) {
        std::string message(
            "Failed to read from the binary file for PostgreSQL "
            "(isTheEnd function in postgres.cpp).");
        throw DataMigratorException(message);
    }
    colNumber = be16toh(colNumber);
    if (colNumber == -1) {
        return true;
    } else if (colNumber == attributes.size()) {
        return false;
    } else {
        std::string message(
            "PostgreSQL input file. Expected number of attributes: "
            + boost::lexical_cast < std::string
            > (attributes.size())
            + "; The number of attributes in the file: "
            + boost::lexical_cast < std::string
            > (colNumber) + ".");
        throw DataMigratorException(message);
    }
}

/* read first 19 bytes */
void Postgres::readFileHeader() {
    size_t headerSize = 19;
    char buffer[headerSize];
    size_t numberOfObjectsRead = fread(buffer, headerSize, 1, fp);
    if (numberOfObjectsRead != 1) {
        std::string message(
            "Failed to read from the binary file for PostgreSQL "
            "(readFileHeader function in postgres.cpp).");
        throw DataMigratorException(message);
    }
    if (strncmp(buffer, BinarySignature, 19) != 0) {
        throw DataMigratorException(
            "Unrecognized PostgreSQL binary signature!");
    }
}

void Postgres::writeFileHeader() {
    // 19 bytes for the signature
    fwrite(&BinarySignature, 19, 1, fp);
}

/* Each tuple begins with a 16-bit integer count of the number of fields in the tuple. */
void Postgres::writeRowHeader() {
    uint16_t colNumberPostgres = htobe16(attributes.size());
    fwrite(&colNumberPostgres, 2, 1, fp);
}

/* The file trailer consists of a 16-bit integer word containing -1.
   This is easily distinguished from a tuple's field-count word. */
void Postgres::writeFileFooter() {
    uint16_t trailer = -1;
    uint16_t trailerPostgres = htobe16(trailer);
    fwrite(&trailerPostgres, 2, 1, fp);
}


text:
//
//  cnn.cpp
//  Jacobian
//
//  Created by David Freifeld
//

#include "bpnn.hpp"
#include "utils.hpp"

//#include <Eigen/unsupported/CXX11/Tensor>
#include "cnn.hpp"

#define LARGE_NUM 1000000 // Remove me.

#if (!RECKLESS)
#define checknan(x, loc) if(x==INFINITY || x==NAN || x == -INFINITY) throw ValueError("Detected NaN in operation", loc)
#else
#define checknan(x, loc)
#endif

// NOTE: Below three functions not mine, from https://compvisionlab.wordpress.com/2014/01/01/c-code-for-reading-mnist-data-set/
int ReverseInt (int i)
{
    unsigned char ch1, ch2, ch3, ch4;
    ch1=i&255;
    ch2=(i>>8)&255;
    ch3=(i>>16)&255;
    ch4=(i>>24)&255;
    return((int)ch1<<24)+((int)ch2<<16)+((int)ch3<<8)+ch4;
}

void ReadMNIST(int NumberOfImages, int DataOfAnImage,std::vector<std::vector<double>> &arr)
{
    arr.resize(NumberOfImages,std::vector<double>(DataOfAnImage));
    std::ifstream file ("./t10k-images-idx3-ubyte",std::ios::binary);
    if (file.is_open())
    {
        int magic_number=0;
        int number_of_images=0;
        int n_rows=0;
        int n_cols=0;
        file.read((char*)&magic_number,sizeof(magic_number));
        magic_number= ReverseInt(magic_number);
        file.read((char*)&number_of_images,sizeof(number_of_images));
        number_of_images= ReverseInt(number_of_images);
        file.read((char*)&n_rows,sizeof(n_rows));
        n_rows= ReverseInt(n_rows);
        file.read((char*)&n_cols,sizeof(n_cols));
        n_cols= ReverseInt(n_cols);
        for(int i=0;i<number_of_images;++i)
        {
            for(int r=0;r<n_rows;++r)
            {
                for(int c=0;c<n_cols;++c)
                {
                    unsigned char temp=0;
                    file.read((char*)&temp,sizeof(temp));
                    arr[i][(n_rows*r)+c]= (double)temp;
                }
            }
        }
    }
}

unsigned char* read_mnist_labels(std::string full_path, int number_of_labels) {
    auto reverseInt = [](int i) {
        unsigned char c1, c2, c3, c4;
        c1 = i & 255, c2 = (i >> 8) & 255, c3 = (i >> 16) & 255, c4 = (i >> 24) & 255;
        return ((int)c1 << 24) + ((int)c2 << 16) + ((int)c3 << 8) + c4;
    };

    typedef unsigned char uchar;

    std::ifstream file(full_path, std::ios::binary);

    if(file.is_open()) {
        int magic_number = 0;
        file.read((char *)&magic_number, sizeof(magic_number));
        magic_number = reverseInt(magic_number);

        if(magic_number != 2049) throw std::runtime_error("Invalid MNIST label file!");

        file.read((char *)&number_of_labels, sizeof(number_of_labels)), number_of_labels = reverseInt(number_of_labels);

        uchar* _dataset = new uchar[number_of_labels];
        for(int i = 0; i < number_of_labels; i++) {
            file.read((char*)&_dataset[i], 1);
        }
        return _dataset;
    } else {
        throw std::runtime_error("Unable to open file `" + full_path + "`!");
    }
}

ConvLayer::ConvLayer(int x, int y, int stride, int kern_x, int kern_y, int pad, std::function<float(float)> activ, std::function<float(float)> activ_deriv)
    :stride_len(stride), padding(pad), activation(activ), activation_deriv(activ_deriv)
{
    pad*=2;
    input = new Eigen::MatrixXf (x+pad,y+pad);
    dZ = new Eigen::MatrixXf (x+pad,y+pad);
    for (int i = 0; i < (x+pad)*(y+pad); i++) {
        (*input)((int)i / (y+pad),i%(y+pad)) = 0;
        (*dZ)((int)i / (y+pad),i%(y+pad)) = 0;
    }
    kernel = new Eigen::MatrixXf (kern_x, kern_y);
    for (int i = 0; i < kern_x*kern_y; i++) {
        (*kernel)((int)i / kern_y,i%kern_y) = (float) rand() / RAND_MAX;
    }
    output = new Eigen::MatrixXf ((x-kern_x+1+pad/stride_len), (y-kern_y+1+pad/stride_len));
    for (int i = 0; i < (x-kern_y+1+pad/stride_len)*(y-kern_x+1+pad/stride_len); i++) {
        (*output)((int)i / (y-kern_y+1+pad/stride_len),i%(y-kern_y+1+pad/stride_len)) = 0;
    }
    bias = 0;
};

void ConvLayer::convolute()
{
    for (int i = 0; i < (*input).rows(); i++) {
        for (int j = 0; j < (*input).cols(); j++) {
            (*dZ)(i, j) = activation_deriv((*input)(i, j));
            (*input)(i, j) = activation((*input)(i, j));
        }
    }
    for (int i = 0; i < output->rows(); i+=stride_len) {
        for (int j = 0; j < output->cols(); j+=stride_len) {
            (*output)(i, j) = (*kernel * (input->block(i, j, kernel->rows(), kernel->cols()))).sum();
        }
    }
    *output = (output->array() + bias).matrix();
}

void ConvLayer::set_input(Eigen::MatrixXf* matrix)
{
    input->block(padding, padding, matrix->rows(), matrix->cols()) = *matrix;
}

// Will eventually be different from ConvLayer
PoolingLayer::PoolingLayer(int x, int y, int stride, int kern_x, int kern_y, int pad)
    :stride_len(stride), padding(pad)
{
    input = new Eigen::MatrixXf (x+pad,y+pad);
    for (int i = 0; i < (x+pad)*(y+pad); i++) {
        (*input)((int)i / (y+pad),i%(y+pad)) = 0;
    }
    kernel = new Eigen::MatrixXf (kern_x, kern_y);
    for (int i = 0; i < kern_x*kern_y; i++) {
        (*kernel)((int)i / kern_y,i%kern_y) = (float) rand()/RAND_MAX;
    }
    output = new Eigen::MatrixXf (x-kern_x+1, y-kern_y+1);
    for (int i = 0; i < (x-kern_x+1)*(y-kern_y+1); i++) {
        (*output)((int)i / (y-kern_y+1),i%(y-kern_y+1)) = (float) rand()/RAND_MAX;
    }
};

void PoolingLayer::pool()
{
  // It doesn't look like anything better than O(n^4) is doable for this as kernel needs to go through matrix and you need to index kernel. LOOK INTO ME!!
    float maxnum = -LARGE_NUM;
    for (int i = 0; i < input->cols() - kernel->cols(); i+=stride_len) {
        for (int j = 0; j < input->rows() - kernel->rows(); j+=stride_len) {
            for (int k = 0; k < kernel->cols(); k++) {
                for (int l = 0; l < kernel->rows(); l++) {
                    if ((input->block(j, i, kernel->rows(), kernel->cols()))(l, k) > maxnum) {
                        maxnum = (input->block(j, i, kernel->rows(), kernel->cols()))(l, k);
                    }
                }
            }
        }
    }
}

ConvNet::ConvNet(const char* path, float learn_rate, float bias_rate, Regularization reg, float l, float ratio)
    :Network(path, 1, learn_rate, bias_rate, reg, l, ratio), preprocess_length{0}
{
    ReadMNIST(10000,784,data);
    data_labels = read_mnist_labels("./t10k-labels-idx1-ubyte",10000);
    labels = new Eigen::MatrixXf (1, 1);
}

void ConvNet::add_conv_layer(int x, int y, int stride, int kern_x, int kern_y, int pad, std::function<float(float)> activ, std::function<float(float)> activ_deriv)
{
    preprocess_length+=1;
    conv_layers.emplace_back(x,y,stride,kern_x, kern_y,pad,activ,activ_deriv);
}

// May make this inaccessible to user code and just have it called from add_conv_layer as pooling is basically always paired with conv.
void ConvNet::add_pool_layer(int x, int y, int stride, int kern_x, int kern_y, int pad)
{
    pool_layers.emplace_back(x,y,stride,kern_x,kern_y,pad);
}

void ConvNet::initialize()
{
    for (int i = 0; i < length-1; i++) {
        layers[i].init_weights(layers[i+1]);
    }
}

void ConvNet::next_batch()
{
    for (int i = 0; i < 784; i++) {
        (*conv_layers[0].input)(i/28, i%28) = data[batches][i];
    }
    (*labels)(0,0) = (float)(int)data_labels[batches];
}

void ConvNet::process()
{
    // Assumes pooling is immediately after any conv layer.
    for (int i = 0; i < preprocess_length-1; i++) {
        conv_layers[i].convolute();
        //   pool_layers[i].input = conv_layers[i].output;
        //   pool_layers[i].pool();
        conv_layers[i+1].input = conv_layers[i].output;
    }
    conv_layers[preprocess_length-1].convolute();
    //pool_layers[preprocess_length-1].input = conv_layers[preprocess_length-1].output;
    //pool_layers[preprocess_length-1].pool();
    //  std::cout << "Output:\n" << *pool_layers[preprocess_length-1].output << "\n\n";
    Eigen::Map<Eigen::RowVectorXf> flattened (conv_layers[preprocess_length-1].output->data(), conv_layers[preprocess_length-1].output->size());
    // std::cout << "Flattened:\n" << flattened << "\n\n";
    for (int i = 0; i < flattened.cols(); i++) {
        (*layers[0].contents)(0, i) = flattened[i];
    }
}

void ConvNet::set_label(Eigen::MatrixXf newlabels)
{
    *labels = newlabels;
}

void ConvNet::list_net()
{
    for (int i = 0; i < preprocess_length; i++) {
        std::cout << "-----------------------\nCONVOLUTIONAL LAYER " << i << "\n-----------------------\n\n\u001b[31mGENERAL INFO:\x1B[0;37m\nStride: " << conv_layers[i].stride_len << "\nPadding: " << conv_layers[i].padding <<  "\n\n\u001b[31mINPUT:\x1B[0;37m\n" << *conv_layers[i].input << "\n\n\u001b[31mKERNEL:\x1B[0;37m\n" << *conv_layers[i].kernel << "\n\n\u001b[31mOUTPUT:\x1B[0;37m\n" << *conv_layers[i].output << "\n\n\u001b[31mBIAS:\x1B[0;37m\n" << conv_layers[i].bias << "\n\n\n";
        //std::cout << "-----------------------\nPOOLING LAYER " << i << "\n-----------------------\n\n\u001b[31mGENERAL INFO:\x1B[0;37m\nStride: " << pool_layers[i].stride_len << "\nPadding: " << conv_layers[i].padding << "\n\n\u001b[31mINPUT:\x1B[0;37m\n" << *pool_layers[i].input << "\n\n\u001b[31mKERNEL:\x1B[0;37m\n-" << *pool_layers[i].kernel << "\n\n\u001b[31mOUTPUT:\x1B[0;37m\n" << *pool_layers[i].output << "\n\n\n";
    }
    std::cout << "-----------------------\nINPUT LAYER (LAYER 0)\n-----------------------\n\n\u001b[31mGENERAL INFO:\x1B[0;37m\nActivation Function: " << layers[0].activation_str << "\n\n\u001b[31mACTIVATIONS:\x1B[0;37m\n" << *layers[0].contents << "\n\n\u001b[31mWEIGHTS:\x1B[0;37m\n" << *layers[0].weights << "\n\n\u001b[31mBIASES:\x1B[0;37m\n" << *layers[0].bias << "\n\n\n";
    for (int i = 1; i < length-1; i++) {
        std::cout << "-----------------------\nLAYER " << i << "\n-----------------------\n\n\u001b[31mGENERAL INFO:\x1B[0;37m\nActivation Function: " << layers[i].activation_str << "\n\n\u001b[31mACTIVATIONS:\x1B[0;37m\n" << *layers[i].contents << "\n\n\u001b[31mBIASES:\x1B[0;37m\n" << *layers[i].bias << "\n\n\u001b[31mWEIGHTS:\x1B[0;37m\n" << *layers[i].weights << "\n\n\n";
    }
    std::cout << "-----------------------\nOUTPUT LAYER (LAYER " << length-1 << ")\n-----------------------\n\n\u001b[31mGENERAL INFO:\x1B[0;37m\nActivation Function: " << layers[length-1].activation_str <<"\n\n\u001b[31mACTIVATIONS:\x1B[0;37m\n" << *layers[length-1].contents << "\n\n\u001b[31mBIASES:\x1B[0;37m\n" << *layers[length-1].bias <<  "\n\n\n";
}

void ConvNet::backpropagate()
{
    list_net();
    char a;
    std::cin >> a;
    std::vector<Eigen::MatrixXf> gradients;
    gradients.push_back(Network::backpropagate());
    Eigen::Map<Eigen::MatrixXf> reshaped(gradients[gradients.size()-1].data(),
                           conv_layers.back().output->rows(),
                           conv_layers.back().output->cols());
    gradients[gradients.size()-1] = reshaped;
    std::vector<Eigen::MatrixXf> conv_deltas;
    for (int layer = conv_layers.size()-1; layer >= 0; layer--) {
        conv_deltas.emplace_back(conv_layers[layer].kernel->rows(),
                                 conv_layers[layer].kernel->cols());
        std::cout << conv_layers[layer].input->cols() << " " << gradients.back().cols() << "\n";
        for (int i = 0; i < conv_layers[layer].input->rows() - gradients.back().rows() + 1; i++) {
            for (int j = 0; j < conv_layers[layer].input->cols() - gradients.back().cols() + 1; j++) {
                conv_deltas[conv_deltas.size()-1](i, j) = (gradients.back() * conv_layers[layer].input->block(i, j, gradients.back().rows(), gradients.back().cols())).sum();
            }
        }
        *conv_layers[layer].kernel -= conv_deltas.back();

        Eigen::MatrixXf flipped_kernel =
            Eigen::MatrixXf::Zero(conv_layers[layer].kernel->rows(), conv_layers[layer].kernel->cols());
        flipped_kernel = conv_layers[layer].kernel->transpose().colwise().reverse().transpose().colwise().reverse();
        Eigen::MatrixXf padded_grad = Eigen::MatrixXf::Zero(gradients.back().rows() + ((flipped_kernel.rows() - 1)*2), gradients.back().cols() + ((flipped_kernel.cols() - 1)*2));
        padded_grad.block(flipped_kernel.rows() - 1, flipped_kernel.cols() - 1, gradients.back().rows(), gradients.back().cols()) = gradients.back();
        Eigen::MatrixXf final_grad (padded_grad.rows() - flipped_kernel.rows() + 1, padded_grad.cols() - flipped_kernel.cols() + 1);
        for (int i = 0; i < padded_grad.rows() - flipped_kernel.rows() + 1; i++) {
            for (int j = 0; j < padded_grad.cols() - gradients.back().cols() + 1; j++) {
                final_grad(i, j) = (flipped_kernel * padded_grad.block(i, j, flipped_kernel.rows(), flipped_kernel.cols())).sum();
            }
        }
        gradients.push_back(final_grad.cwiseProduct(*conv_layers[layer].dZ));
    }
    list_net();
    assert(2<1);
}

void ConvNet::train()
{
    float cost_sum = 0;
    float acc_sum = 0;
    for (int i = 0; i <= 100; i++) {
        if (i != instances-batch_size) { // Don't try to advance batch on final batch.
            next_batch();
        }
        process();
        feedforward();
        backpropagate();
        cost_sum += cost();
        acc_sum += accuracy();
        batches++;
    }
    list_net();
    epoch_acc = 1.0/(100) * acc_sum;
    epoch_cost = 1.0/(100) * cost_sum;
    printf("Epoch %i complete - cost %f - acc %f\n", epochs, epoch_cost, epoch_acc);
    batches=0;
    decay();
    epochs++;
}

int main()
{
    ConvNet net ("../data_banknote_authentication.txt", 0.05, 0.01, L2, 0, 0.9);
    Eigen::MatrixXf labels (1,1);
    net.add_conv_layer(28, 28, 1, 9, 9, 0, lecun_tanh, lecun_tanh_deriv);
    //  net.add_pool_layer(20,20,1,6,6,0);
    net.add_conv_layer(20, 20, 1, 6, 6, 0, lecun_tanh, lecun_tanh_deriv);
    std::cout << net.conv_layers[net.conv_layers.size()-1].output->rows() << "\n";
    //net.add_pool_layer(10,10,1,2,2,0);
    net.add_layer(400, "sigmoid", sigmoid, sigmoid_deriv);
    net.add_layer(5, "lecun_tanh", lecun_tanh, lecun_tanh_deriv);
    net.add_layer(10, "resig", rectifier(sigmoid), rectifier(sigmoid_deriv));
    //  net.list_net();
    //  net.init_decay("step", 1, 2);
    net.initialize();
    //net.list_net();

    for (int i = 0; i < 1; i++) {
        net.train();
    }
    net.list_net();
}

max_stars_repo_name: richardfeynmanrocks/ml-in-parallel

max_issues_repo_name: quantumish/Jacobian

max_forks_repo_name: quantumish/Jacobian


text:
#include "leveldb_storage.h"

#include <boost/chrono/system_clocks.hpp>
#include <boost/filesystem/operations.hpp>
#include <boost/throw_exception.hpp>

#include <leveldb/comparator.h>

namespace memhook {
  namespace {
    class SystemClockComparator : public leveldb::Comparator {
    public:
      int Compare(const leveldb::Slice &lhs, const leveldb::Slice &rhs) const {
        const chrono::system_clock::rep lhs_rep =
                *reinterpret_cast<const chrono::system_clock::rep *>(lhs.data());

        const chrono::system_clock::rep rhs_rep =
                *reinterpret_cast<const chrono::system_clock::rep *>(rhs.data());

        if (lhs_rep < rhs_rep)
          return -1;

        if (lhs_rep > rhs_rep)
          return 1;

        return 0;
      }

      const char *Name() const {
        return "SystemClockComparator";
      }
      void FindShortestSeparator(std::string *, const leveldb::Slice &) const {}
      void FindShortSuccessor(std::string *) const {}
    };

    SystemClockComparator SYSTEM_CLOCK_COMPARATOR;
  }  // namespace

  LevelDBStorage::LevelDBStorage(const char *path, std::size_t cache_size_mb)
      : m_path(path)
      , m_cache_size_mb(cache_size_mb) {
    if (!boost::filesystem::is_directory(m_path)) {
      if (boost::filesystem::exists(m_path))
        boost::filesystem::remove(m_path);
      boost::filesystem::create_directory(m_path);
    }

    leveldb::Options options;
    options.create_if_missing = true;
    options.compression = leveldb::kNoCompression;

    m_symtab = OpenDB(options, m_path / "symtab");
    m_shltab = OpenDB(options, m_path / "shltab");

    options.comparator = &SYSTEM_CLOCK_COMPARATOR;
    if (m_cache_size_mb) {
      options.block_cache = leveldb::NewLRUCache(m_cache_size_mb * 1024 * 1024);
      m_dcache = unique_ptr<leveldb::Cache>(options.block_cache);
    }
    m_dattab = OpenDB(options, m_path / "dattab");
  }

  unique_ptr<leveldb::DB> LevelDBStorage::OpenDB(const leveldb::Options &options,
          const boost::filesystem::path &path) const {
    leveldb::DB *db = NULL;
    leveldb::Status status = leveldb::DB::Open(options, path.string(), &db);
    if (!status.ok()) {
      BOOST_THROW_EXCEPTION(std::runtime_error(status.ToString()));
    }
    return unique_ptr<leveldb::DB>(db);
  }

  void LevelDBStorage::Add(uintptr_t address,
          std::size_t memsize,
          const CallStackInfo &callstack,
          const chrono::system_clock::time_point &timestamp) {}

  bool LevelDBStorage::Remove(uintptr_t address) {
    return false;
  }

  bool LevelDBStorage::UpdateSize(uintptr_t address, std::size_t memsize) {
    return false;
  }

  void LevelDBStorage::Clear() {}

  void LevelDBStorage::Flush() {}

  std::string LevelDBStorage::GetName() const {
    return m_path.string();
  }

  unique_ptr<MappedStorage> NewLevelDBStorage(const char *path, std::size_t cache_size_mb) {
    return unique_ptr<MappedStorage>(new LevelDBStorage(path, cache_size_mb));
  }

}  // memhook

max_stars_repo_name: sergzub/memhook

max_issues_repo_name: sergzub/memhook

max_forks_repo_name: sergzub/memhook


text:
//  (C) Copyright 2014 Alvaro J. Genial (http://alva.ro)
//  Use, modification and distribution are subject to the Boost Software License, Version 1.0.
//  (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt).

#ifndef AJG_SYNTH_ADAPTERS_BASE_ADAPTER_HPP_INCLUDED
#define AJG_SYNTH_ADAPTERS_BASE_ADAPTER_HPP_INCLUDED

#include <ajg/synth/support.hpp>

#include <string>
#include <istream>
#include <ostream>
#include <typeinfo>
#include <stdexcept>
#include <functional>

#include <boost/foreach.hpp>
#include <boost/optional.hpp>

#include <ajg/synth/exceptions.hpp>

namespace ajg {
namespace synth {

using boost::optional;

template <class Behavior, class Adapted>
struct adapter;

//
// base_adapter
////////////////////////////////////////////////////////////////////////////////////////////////////

template <class Behavior>
struct base_adapter {
  public:

    typedef Behavior                                                            behavior_type;
    typedef base_adapter                                                        base_type;
    typedef typename behavior_type::traits_type                                 traits_type;
    typedef typename behavior_type::value_type                                  value_type;
    typedef typename behavior_type::adapter_type                                adapter_type;

    typedef typename traits_type::char_type                                     char_type;
    typedef typename traits_type::size_type                                     size_type;
    typedef typename traits_type::boolean_type                                  boolean_type;
    typedef typename traits_type::floating_type                                 floating_type;
    typedef typename traits_type::string_type                                   string_type;
    typedef typename traits_type::datetime_type                                 datetime_type;
    typedef typename traits_type::duration_type                                 duration_type;
    typedef typename traits_type::istream_type                                  istream_type;
    typedef typename traits_type::ostream_type                                  ostream_type;

    typedef typename value_type::iterator                                       iterator;
    typedef typename value_type::const_iterator                                 const_iterator;
    typedef typename value_type::range_type                                     range_type;

  public:

    inline base_adapter() {}
    virtual ~base_adapter() {}

  public:

    // TODO: Instead of throwing invalid_method here, defer default behavior to Traits.

    virtual std::type_info const& type() const = 0;

    virtual boolean_type  is_numeric()  const { return false; }

    virtual floating_type to_floating()   const { AJG_SYNTH_THROW(invalid_method("to_floating")); }
    virtual boolean_type  to_boolean()  const { AJG_SYNTH_THROW(invalid_method("to_boolean")); }
    virtual datetime_type to_datetime() const { AJG_SYNTH_THROW(invalid_method("to_datetime")); }
    virtual string_type   to_string()   const { return behavior_type::to_string(*this); }

    virtual void input (istream_type& in)        { AJG_SYNTH_THROW(invalid_method("input")); }
    virtual void output(ostream_type& out) const { AJG_SYNTH_THROW(invalid_method("output")); }

    // TODO: Refactor these into a single to_range function.
    virtual const_iterator begin() const { AJG_SYNTH_THROW(invalid_method("begin")); }
    virtual const_iterator end()   const { AJG_SYNTH_THROW(invalid_method("end")); }

    virtual const_iterator       find(value_type const& value) const { AJG_SYNTH_THROW(invalid_method("find")); }
    virtual optional<value_type> index(value_type const& what) const { AJG_SYNTH_THROW(invalid_method("index")); }

  private:

    template <class T, class A>                    friend struct synth::adapter;
    template <class B, class T, class A, class S>  friend struct forwarding_adapter;
    template <class T, template <class> class V>   friend struct value_facade;

    template <class Adapter> // TODO: Deal with forwarding_adapters.
    inline Adapter const* get() const { return dynamic_cast<Adapter const*>(this); }

    template <class T> // TODO: Deal with forwarding_adapters.
    inline T const& get_adapted() const {
        typedef synth::adapter<Behavior, T> specialized_type;
        specialized_type const* const specialization = this->template get<specialized_type>();
        BOOST_ASSERT(specialization);
        return specialization->adapted_;
    }

  protected:

    virtual boolean_type equal_adapted(adapter_type const& that) const = 0;
    virtual boolean_type less_adapted (adapter_type const& that) const = 0;

    template <class Adapter>
    inline boolean_type equal_as(adapter_type const& that) const {
        Adapter const* const this_ = this->template get<Adapter>();
        Adapter const* const that_ = that->template get<Adapter>();
        return this_ != 0 && that_ != 0 && std::equal_to<typename Adapter::bare_adapted_type>()(this_->adapted_, that_->adapted_);
    }

    template <class Adapter>
    inline boolean_type less_as(adapter_type const& that) const {
        Adapter const* const this_ = this->template get<Adapter>();
        Adapter const* const that_ = that->template get<Adapter>();
        return this_ != 0 && that_ != 0 && std::less<typename Adapter::bare_adapted_type>()(this_->adapted_, that_->adapted_);
    }
};

}} // namespace ajg::synth

#endif // AJG_SYNTH_ADAPTERS_BASE_ADAPTER_HPP_INCLUDED

max_stars_repo_name: legutierr/synth

max_issues_repo_name: legutierr/synth

max_forks_repo_name: legutierr/synth


text:
/*
@copyright Louis Dionne 2014
Distributed under the Boost Software License, Version 1.0.
(See accompanying file LICENSE.md or copy at http://boost.org/LICENSE_1_0.txt)
 */

#include <boost/hana/list/mcd.hpp>

#include <boost/hana/detail/assert.hpp>
#include <boost/hana/detail/constexpr.hpp>
#include <boost/hana/detail/minimal/comparable.hpp>
#include <boost/hana/detail/minimal/list.hpp>
using namespace boost::hana;

template <int i>
constexpr auto x = detail::minimal::comparable<>(i);

template <typename mcd>
void test() {
    BOOST_HANA_CONSTEXPR_LAMBDA auto list = detail::minimal::list<mcd>;
    auto p = list; // tests are unreadable otherwise. mnemonic: pair

    BOOST_HANA_CONSTANT_ASSERT(unzip(list(p())) == list());

    BOOST_HANA_CONSTANT_ASSERT(unzip(list(p(), p())) == list());
    BOOST_HANA_CONSTEXPR_ASSERT(unzip(list(p(x<0>, x<2>), p(x<1>, x<3>, x<4>))) == list(p(x<0>, x<1>), p(x<2>, x<3>)));

    BOOST_HANA_CONSTANT_ASSERT(unzip(list(p(),     p(),     p()))     == list());
    BOOST_HANA_CONSTANT_ASSERT(unzip(list(p(x<0>), p(),     p()))     == list());
    BOOST_HANA_CONSTANT_ASSERT(unzip(list(p(),     p(x<1>), p()))     == list());
    BOOST_HANA_CONSTANT_ASSERT(unzip(list(p(),     p(),     p(x<2>))) == list());
    BOOST_HANA_CONSTANT_ASSERT(unzip(list(p(x<0>), p(x<1>), p()))     == list());
    BOOST_HANA_CONSTANT_ASSERT(unzip(list(p(),     p(x<1>), p(x<2>))) == list());
    BOOST_HANA_CONSTANT_ASSERT(unzip(list(p(x<0>), p(),     p(x<2>))) == list());
    BOOST_HANA_CONSTEXPR_ASSERT(unzip(list(p(x<0>), p(x<1>), p(x<2>))) == list(p(x<0>, x<1>, x<2>)));

    BOOST_HANA_CONSTEXPR_ASSERT(unzip(list(p(x<0>, x<3>), p(x<1>, x<4>), p(x<2>, x<5>))) == list(p(x<0>, x<1>, x<2>), p(x<3>, x<4>, x<5>)));
}

int main() {
    test<List::mcd<void>>();
}

max_stars_repo_name: rbock/hana

max_issues_repo_name: rbock/hana

max_forks_repo_name: rbock/hana